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Summary 

This work is concerned with modelling the mechanical deformation or constitutive behavior 
of composites comprised of a periodic microstructure under small displacement conditions at 
elevated temperature. A mesomechanics approach [1] is adopted which relates the microme- 
chanical behavior of the heterogeneous composite with its in-service macroscopic behavior. 

Two different methods, one based on a Fourier series approach and the other on a Green’s 
function approach, are used in modelling the micromechanical behavior of the composite 
material. Although the constitutive formulations are based on a micromechanical approach, 
it should be stressed that the resulting equations are volume averaged to produce overall 
“effective” constitutive relations which relate the bulk, volume averaged, stress increment to 
the bulk, volume averaged, strain increment. As such, they are macromodels which can be 
used directly in nonlinear finite element programs such as MARC, ANSYS and ABAQUS or 
in boundary element programs such as BEST3D. 

In developing the volume averaged or “effective” macromodels from the micromechanical 
models, both approaches (i.e. Fourier series and Green’s function) will require the evalua- 
tion of volume integrals containing the spatially varying strain distributions throughout the 
composite material. By assuming that the strain distributions are spatially constant within 
each constituent phase — or within a given subvolume within each constituent phase — of the 
composite material, the volume integrals can be obtained in closed form. This simplified 
micromodel can then be volume averaged to obtain an “effective” macromodel suitable for 
use in the MARC, ANSYS and ABAQUS nonlinear finite element programs via user consti- 
tutive subroutines such as HYPELA and CMUSER. This “effective” macromodel can be used 
in a nonlinear finite element structural analysis to obtain the strain-temperature history at 
those points in the structure where thermomechanical cracking and damage are expected to 
occur, the so called “damage critical” points of the structure. The “exact” micromechanical 
models can then be subjected to the overall “effective” strain-temperature history obtained 
at the “damage critical” location and used outside of the finite element program to evaluate 
the heterogeneous stress-strain history throughout each constituent phase of the composite 
material. This variation must be known in order to evaluate the damage history variation 
throughout each constituent phase of the composite material. 

‘Work funded by NASA Grant NAG.3-882. 
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1 Introduction 


The ultimate objective of this work is to produce a computer program to analyze the hetero- 
geneous stress and strain history variation at the “damage critical” locations of a composite 
structure operating at elevated temperature. This report describes some of the theoreti- 
cal foundations for the program. A mesomechanics [1] approach is adopted which relates 
the micromechanical behavior of the heterogeneous composite to its in-service macroscopic 
behavior. 

Some composites are actually comprised of a periodic microstructure whilst others are 
possessed of an essentially randomly distributed microstructure. Pictures of metal matrix 
composites (tungsten-fiber-reinforced superalloys) which exhibit a periodic microstructure 
are shown in Fig. 1 which is taken from the article by Petrasek et al [2]. When the fibers in a 
composite material occupy a large volume fraction of the material, the induced deformation 
in one fiber interacts with and alters the induced deformation in the neighboring fibers. When 
the fibers are densely packed the interaction effect becomes dominant and must be accounted 
for in the constitutive formulation. 

At NASA-Lewis Chamis and his colleagues [3,4,5] employ two different approaches for 
analyzing the behavior of structural composites. One method employs a sophisticated finite 
element analysis of a periodic microstructure. A unit cell in the periodic microstructure is 
modelled with one hundred and ninety two three-dimensional elements and the eight nearest 
neighbor cells of the fibrous composite are modelled with superelements. By applying the 
strain-temperature history at the “damage critical” location in the composite structure to 
the superelement model, the stress-strain history throughout the unit cell can be computed 
and used to estimate the maximum damage in the composite structure. This method will 
necessarily require large resources in computer time and memory to analyze the viscoplastic 
behavior of the composite structure under in-service thermomechanical loading conditions. 

Another approach adopted by Chamis and his colleagues [3,4,5] — which avoids large com- 
puter resources — is to employ composite micromechanics theory to derive simplified rela- 
tionships which describe the thermomechanical constitutive behavior of multilayered fibrous 
composites. 

When suitable boundary conditions are applied to the superelement model of the periodic 
microstructure, it is possible to predict the elastic properties of the equivalent homogenized 
material. A comparison [5] with the homogenized elastic properties predicted by the simplified 
micromechanical equations generally shows good agreement with the superelement model 
except for the Poisson ratios. At high volume fractions (~ 60%) the longitudinal Poisson’s 
ratio for unidirectional fibers predicted by the simplified equations is about 15% too small, 
whilst the transverse Poisson’s ratio is about 30% too small. These anomalies occur because 
the interaction between the fibers is not accounted for in the simplified micromechanical 
model. This may be important when considering the highly nonlinear behavior of viscoplastic 
composites at elevated temperature. 

Dvorak [6] and Dvorak and Bahei-El-Din et al [7,8,9] have also made great progress in 
modelling the micromechanical behavior of nonlinear composite materials and are embarked 
on a combined experimental and theoretical effort. The variation of the stress-strain history 
throughout the unit cell of a periodic microstructure is obtained with a finite element analysis 
in which the interaction effects of the surrounding cells is accounted for by applying periodic 
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boundary conditions to the surface of the unit cell. 

Work on the theoretical foundations behind the homogenization of micromechanical con- 
stitutive models to produce bulk macroscopic models has been under way in France by Devries 
and Lene [10], Lene [11,12], Duvaut [13], Renard and Marmonier [14], Lene and Leguillon 
[15] and Sanchez-Palencia [16,17]. These references give a good account of the work be- 
ing conducted in France by other researchers. Much of this work is founded on the use of 
multivariable asymptotic techniques [18]. In an infinite periodic structure the stress-strain 
history in each unit periodic cell is, perforce, identical. Due to the finite size of the compos- 
ite structure the effects of surface tractions and displacements on the surface will cause the 
stress-strain history to vary from cell to cell. If the unit cell is much smaller than the size 
of the structure this variation from cell to cell will be small. If L is a typical dimension of 
the unit periodic cell and D is a typical dimension of the composite structure, then the ratio 
L/D is a small parameter of the problem. The displacement variation throughout the unit 
cell will then depend on six spatial variables, i.e., 

U{ Tl'i i^X i , X2 , X$ 5 ^3) 

where x * = XtL/D. The spatial variables x* take into account the slow variation of the 
displacement from cell to cell due to the finite size of the ratio L/D when it* is a periodic 
function of the variables x,. By expanding the displacement and other spatial variables of the 
problem into a series in powers of L/D and equating like powers in the perturbation expansion, 
it is possible to obtain the effect of the finite size of the structure on the deformation behavior 
in the unit cell. Due to the perturbative assumption of small L/D this method is not expected 
to be valid for thin composite sections or to be applicable at those places in the structure 
where surface effects or nonperiodic inclusions are important. 

Rather than employing finite element techniques to determine the stress-strain history 
variation throughout the unit periodic cell, Aboudi [19] has recently developed a macro- 
scopic formulation for periodic composites based on volume averaging a viscoplastic consti- 
tutive model over the unit periodic cell. This work expands the heterogeneous displacement 
throughout the constituent phases of the unit cell as linear and higher order functions of the 
coordinates. Good agreement with experimental results was achieved by volume averaging 
Bodner’s [20] viscoplastic constitutive model over the unit periodic cell, but the method is 
general and any constitutive model may be used to represent the deformation behavior of the 
constituent phases. The limitation here is that large spatial gradients in the strain history 
may not be accurately modelled by linear or quadratic interpolation functions on the unit 
periodic cell. 

Weng and his colleagues [21,22] have employed self-consistent methods to study the effect 
of inclusion size and volume fraction on the stress distribution in and around spheroidal 
inclusions embedded in an “effective” non-uniform matrix material, and the effect which this 
has on the overall “effective” macroscopic constitutive behavior of the composite. In the 
first paper they point out that the derivation of the fictitious body forces which represent 
the inelastic behavior of the heterogeneous composite material should be obtained from first 
principles rather than using their heuristic approach. In the second paper the Mori- Tanaka 
theorem [24] is used to represent the effect of the heterogeneous composite, and a similar 
procedure is followed in the present work to develop a self-consistent method for composites 
which exhibit a periodic microstructure. In addition, the present report also derives the 


2 



fictitious body forces for a periodic microstructure from first principles. In reference [23] 
Zhu and Weng have used a combined micromechanics and continuum theory approach to 
develop a creep deformation model for particle-strengthened metal matrix composites. They 
stress the fact that the creep resistance of the composite is underestimated when simplified 
metallurgical and mechanics approaches are adopted. 

A comprehensive application of micromechanics to mechanical deformation problems is 
given by Mura [24] in his book “Micromechanics of Defects in Solids” . This work was used by 
Nemat-Nasser and his colleagues who have exploited the mathematical simplicity of a periodic 
microstructure in order to develop elastic, plastic and creep constitutive models [25,26,27,28] 
for composite materials. The assumption of periodicity allows the heterogeneous stress, strain 
and displacement fields to be expanded in a Fourier series, which greatly simplifies the ensuing 
computations. This technique fully accounts for the interaction effects between neighboring 
fibers. Even when the composite is comprised of closely packed fibers distributed at random 
the method gives accurate results [25] for the “effective” elasticity tensor. When densely 
packed fibers form a large volume fraction of the composite material these interaction effects 
play a dominant role and must be included in the calculations. It appears that inclusion of the 
interaction effects can be as, or more, important than inclusion of the random nature of the 
microstructure when the fibers occupy a large volume fraction of the composite material. In 
this report we have developed the Fourier series approach in order to handle the viscoplastic 
behavior of the constituents in the unit periodic cell. 

The nonlinear constitutive behavior of composites with a periodic microstructure can also 
be treated with a Green’s function approach [29,30,31,32,33]. Here, the periodic heteroge- 
neous material property variation — due to the fibers — is treated as a fictitious body force in 
the matrix material. The Green’s function is used to evaluate the displacement due to a unit 
point force in the matrix material and the actual displacement at any point in the composite 
can then be determined by summing (integrating) the effect due to a volume distribution 
of fictitious periodic body forces. It is shown in Appendix B that this method is exactly 
equivalent to the Fourier series approach by invoking a mathematical technique known as the 
Poisson sum formula. The Green’s function approach is more general in that the method can 
also handle the nonperiodic case where there may be inclusions in one unit cell but not in 
the neighboring cells. It is also able to handle surface effects, although the surface integrals 
which represent the surface effects in the Green’s function method could be expanded in a 
Fourier series for thin composite sections. 

The approach adopted in the present work is to develop homogenization techniques which 
can provide simplified macromodels for use in a nonlinear finite element program, similar in 
spirit to the simplified models used at NASA-Lewis, but which account for the viscoplastic 
interaction effects in the periodic structure and which allow surface effects for thin struc- 
tures to be taken into account. Once the strain- temperature history at the “damage critical” 
location has been found from the finite element analysis, it can be used to “drive” the mi- 
cromechanical relations in order to obtain the stress-strain history variation throughout the 
unit cell. These micromechanical relations are the same relations which are used to obtain 
the simplified homogenized constitutive model. When the unit cell is chosen to have the form 
shown in Fig. 2, it is clear that a periodic arrangement of such a microstructure allows for 
the analysis of laminated composite structures. 
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2 Overview of Theoretical Modelling Approaches 

2.1 Outline of Approach 

In order to develop a homogeneous macroscopic constitutive model from micromechanical 
principles it is necessary to know the stress-strain history throughout the unit cell of a periodic 
composite. Some of the approaches which are presently being given currency are described 
in the introduction. 

The Fourier series and Green’s function approaches can be used to compute the viscoplas- 
tic stress-strain history throughout the unit cell of a periodic composite and can be simplified 
to produce a model suitable for use in a nonlinear finite element program. In this report 
the Fourier series and Green’s function approaches are developed and shown to be equivalent 
to each other by means of the Poisson sum formula. This equivalence holds only for an in- 
finitely extended medium. When the medium has a finite size the effect of spatially varying 
displacements and tractions on the surface of the medium must be accounted for. This is 
easily accomplished with the Green’s function method by retaining the appropriate surface 
integral contributions which are discarded in the case of an infinite medium. In the Fourier 
series approach the surface integral could be included, and, in the case of a thin composite 
section which has an infinite surface the integral can be expanded into a Fourier series. In 
fact, the methods can be combined, so that if inclusions are present in one unit cell and not 
in the neighboring cells, their effect can be taken into account in the Fourier series method 
by treating them with a Fourier integral or Green’s function approach. 

An overview of the present work is depicted in Fig. 3. Simplified versions of the mi- 
cromechanical constitutive equations can be volume averaged to produce a macroscopic ho- 
mogenized viscoplastic model. This can then be used in a nonlinear finite element program 
to analyze the structural behavior of a composite structure under in-service thermomechan- 
ical loading conditions. The finite element analysis yields the strain-temperature history at 
the “damage critical” location and this history can then be applied to the micromechanical 
equations to determine the heterogeneous stress-strain history throughout the unit cell. 

A detailed flow chart in Fig. 4 shows the anticipated structural analysis procedure. Both 
the Fourier series and Green’s function approaches can be used to create a coarse subvolume 
model. This coarse model can then be homogenized and included in a user defined constitutive 
subroutine in a nonlinear finite element program. The Green’s function approach can also be 
used to derive a simple self-consistent model for use in the user subroutine. 

2.2 Homogenized Macroscopic Equations 

A periodic composite material is supposed acted upon by an imposed strain increment Ae°- 
and responds in bulk with a stress increment Act®-. These values are then equated to the 
respective volume averaged quantities in order to obtain the “effective” constitutive relation 
for the composite material, i.e., 

A <4 = fJJ Acr b( r ) dV ( r ) and A 4 = \ JJJ Ae p( r ) dV ( r ) (!) 

' v v 

where V is the volume of the body. 
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In section 4 it is shown that the volume averaged or “effective” constitutive relation for 
the composite material can be written as 

H -yJJJ { M - 6D ljU (r) [A^ (r) - A c„ (r)] } dV( r) (2) 

C V c 

where V c is the volume of a unit periodic cell in the composite material, Ae kl (r) is the total 
strain increment at point r in the periodic cell due to the imposed uniform total strain 
increment Ae kl at the surface of the composite, and Ac*/(r) is the strain increment at point 
r in the periodic cell representing the deviation from isothermal elastic behavior. The fourth 
rank tensor r) is defined by the relation 

r) = >?(r) (d' jm - D% u ) (3) 

where i?(r) = 1 in the fiber and $(r) = 0 in the matrix, with D-j k[ denoting the elasticity 
tensor of the fiber and D™ kl that of the matrix. 

In the expression for the average or “effective” constitutive relation in equation 2, the 
quantities Ae° z , D™ kl and 6Dij kl (r ) are given. The deviation strain increment Acjy(r) can 
be obtained throughout the periodic cell as a function of position r by using an explicit 
forward difference method since the stress and state variables in a viscoplastic formulation 
will be known functions of position at the beginning of the increment. Everything is therefore 
known explicitly except the total strain increment Ae kl (r). 


2.3 Fourier Equation Overview 

In the Fourier series approach described in section 4 we find that the total strain increment 
is determined by solving the integral equation, 

i dhoo 

te T kl (r) = te° kl + EE' 9 kiij (C) x 

V c n p = 0 

x /// e ‘ £ (r " r,> { D "rACr. (O - 6D ijr , (r') [A £ (r') - Ac r » (r')]} d.V( r') (4) 

v c 


where the fourth rank tensor g k nj (£) is given by 


««« (0 = \ «) + GW ( O ) 


in which the Christoffel stiffness tensor (£), with inverse M y 1 (£) is defined (c/. [33]) by 
the relation, 

Mij (C) = D™ q C P C„ ( 6 ) 

with ( p = £ p / \/£ m £ m = £ p /£ being a unit vector in the direction of the Fourier wave vector £, 
and £ = denoting the magnitude of the vector £. In equation 4 the sum is taken over 

integer values in which 

27 mi 27 rn 2 . 27 rn 3 . 

SI — — p » S2 7 , S3 " —j {<) 

Li\ L 2 1^3 
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and L\, L 2 , L 3 are the dimensions of the unit periodic cell in the X\, x 2 , X3 directions, so that 
V r = L\ L2L3. The values of n b n 2 , n 3 are given by 

n p = 0, ±1, ±2, ±3, . . . , etc., for p= 1,2,3 (8) 

and the prime on the triple summation signs indicates that the term with ni = n 2 = n$ = 0 
is excluded from the sum. 


2.4 Green’s Equation Overview 

In the Green’s function approach the total strain increment Ae^(r) is determined by solving 
a different integral equation, viz., 


Ae£(r) = A 4 + JJJ U Umn (r - r') {zC„ r .Ac„ (r') - 

V 

- SDmnrs (0 [A# M ~ Ac,, (r')] } dV( t') (9) 


where the fourth rank tensor Ukimn (r — r') gives the kl component of the total strain incre- 
ment at point r due to the ran component of a stress increment applied at point r' in the 
infinite matrix with elasticity tensor £>™ nrs , i.e., 


U^lmn (r 


1 / (PGkm (r - r') + &Gim (r-rQ \ 

2 \ dxidx n dx k dx n J 


( 10 ) 


and the volume integration in equation 9 extends over all the periodic cells in the composite 
material, i.e.. over the entire composite. 

The Green’s function tensor is defined in Appendix A, equation A. 26, by the Fourier 
integral [24,32,33] 


OO 

G„ (r - r') = HI 

— oo 


d 3 K MjJ 1 (C) 
(2tt) 3 K 2 


e -i K.(r-r') 


( 11 ) 


in which the tensor £ is now defined by the relation Q = Ki/K with K = y 'K q K q denoting 
the magnitude of the vector K = (K\, K 2 , K3). 

In Appendix B it is shown, by applying the Poisson sum formula, that equations 4 and 9 
are identical, although the summation extends over the integer values ri\, n 2 , in equation 4 
and extends over the periodic cells in equation 9. 


2.5 Integration of the Equations 

Both equations 4 and 9 are implicit integral equations for the determination of the total strain 
increment Ae^(r), since this unknown quantity appears both on the left hand sides of the 
equations and on the right hand sides under the volume integrations. 

The “effective” constitutive relation given in equation 2 and the total strain increment 
relation, given by either equation 4 or 9, contain the volume integration of the deviation 
strain increment Ac/,.j(r). In the periodic cell the deviation strain increment at any point r 
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will be determined from a unified viscoplastic constitutive relation [34] appropriate to the 
constituent phase in which the point r resides. If a constituent phase is included at the fiber- 
matrix interface, a constitutive relation can also be proposed for this phase, and the resulting 
inelastic strain increment determined for inclusion in the volume integrals. This may be 
important for metal matrix composites where there can be chemical reactions between fiber 
and matrix at elevated temperatures, and for composites where the fiber has been coated to 
enhance overall composite properties. 

Equations 2, 4 and 9 form the basic incremental constitutive equations for determining the 
“effective” overall deformation behavior of a composite material with a periodic microstruc- 
ture. In order to update the stress state in each of the constituent phases in preparation 
for integrating the “effective” constitutive relation over the next increment, the constitutive 
relation, 

Aa^r) = D ijk i( r) (Aej[j(r) - Ac fc ,(r)) (12) 

is used, where r) = D{j kl or D™ kl according as the point r is in the fiber or matrix. 
This relation is used to update the stress cr, ; (r) and, in turn, the internal viscoplastic state 
variables <&( r) at each point r in preparation for computing Ac*;/(r) in the next increment. 

The derivation of the preceding equations and some methods for their solution are dis- 
cussed in the succeeding sections of this report. Numerical solutions will be obtained during 
the research effort from appropriate FORTRAN computer programs. 


3 Periodic Microstructure 


3.1 Volume Averaging 


The periodic composite is supposed acted upon at its surface by a spatially linear displacement 
increment, A«f(r), given by 

Au°(r) = Xj Ae% + XjAuij (13) 

where A and A ufj are the spatially uniform strain and rotation increments at the surface 
of the composite. 

If the matrix material was homogeneous and had no fibers embedded in it, the strain 
increment would be homogeneous and given by 


A4 = i 

1 ] 2 


a(A«?) | 

dxj dxi I 


(14) 


Since this is constant, we may trivially volume average Ae- over the volume V of the homo- 
geneous matrix material to obtain 


A £ ° = I fffl ( a ( AM ' (r) ) + a ( A "°< r >) ) 
« VJJJ 2 l Bxj dx, j 

which, by Gauss’ divergence theorem, may be written as 

A 2 («j( r )Au?( r ) + ni(r)Au°(r)) 

.S’ 


dV(r) 


dS( r) 


(15) 


(16) 
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where the integral extends over the surface of the material and n,(r) denotes the outwardly 
directed unit normal vector at point r on the surface. Thus, by applying the displacement 
increment A«°(r) in equation 13 over the surface of the material to produce the surface strain 
increment given in equation 16, equations 15 and 16 show that the strain increment in the 
matrix material is spatially uniform. 

If the displacement increment A«-(r) in equation 13 is applied to the actual composite 
material, the total displacement increment within the material, A uf (r), will vary in a periodic 
manner due to the assumed geometric periodicity of the composite material, so that 

A uj (r) = Auf(r) + Au*(r) (17) 

where Au°(r) is the displacement increment which would be induced in the homogeneous 
matrix if the fiber phase were absent, and Atq(r) is the perturbation or deviation from the 
homogeneous value due to the presence of the fibers. 

Corresponding to these displacement increments, the total strain increment at any point 
r in the composite, Ae^(r), is given by the relation 


A£ h( r ) = A 4i + Ae*/( r) 


(18) 


where 


A 4t = x 


d(Airg) d(A 4Y] 
dxi dx k ) 


and 


&£ki(r) = - 


( djAuk ) <9(Au,) \ 

\ dxi dx k ) 


(19) 


with Aeh representing the spatially constant total strain increment which would be produced 
on the surface and in the interior of the homogeneous matrix if the fibers were absent, and 
Ae*./(r) representing the deviation from the uniform value due to the presence of the fibers. 
Both the total strain increment Ae^(r) and the perturbed strain increment Ae^(r) vary 
throughout the composite in a periodic manner. 

We define the volume averaged stress and strain increments as Act? and Ae° , respectively. 
The required “effective” constitutive equation for the composite material is then an expression 
relating the volume averaged stress and strain increments. For a function /( r) which varies 
with position the volume average is defined by the relation 

(f) = ^JJ[md v (r) (20) 

v 

Since the composite is assumed to be comprised of a periodic aggregate of identical unit cells, 
we may write 

{f ) = vJJJ f{r)dV (r) (21) 

c V c 

where V c denotes the volume of the unit periodic cell. 

If we volume average the total strain increment in equation 18, we obtain 

<Ae£) = yfff A£ L( r ) dV ( r ) = A 4i + yJJJ A£ w ( r ) dV ( r ) ( 22 ) 

c v r c v c 
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or 


(A 4) = AtJ, + (A<=„> (23) 

But the volume averaged total strain increment is defined as Ash , so that (Ae' kl ) = Ae kl and 

<Ae«) = 0 (24) 

which shows that the volume averaged perturbation strain increment, Ae^r), is equal to 
zero. 


3.2 Eigenstrain and Deviation Strain Increments 

If the elasticity tensor is denoted by Dij k i( r) and the inelastic strain tensor by e kl ( r), then 
the constitutive equation at any point r in the composite material can be written as 

a ij ( r ) = Dijkfir) (e«(r) - eQ(r) - a M (r) (T - T 0 )) (25) 

where a k i( r) is the coefficient of thermal expansion. 

The incremental form of Hooke’s law is 

Affy(r) = Dijki(r) (Ae£(r) - Ac«(r)) (26) 

where Ac k i(r ) denotes the incremental strain representing the deviation from isothermal elas- 
tic conditions and is given by 

Ac w (r) = Aejy(r) + a«(r)Ar (27) 


in which 


Q w( r )A7 = a w (r)AT + (r-T 0 )Aa w (r)- 

- D ki!j( r ) AD ijmn (r) (e£ n ( r ) - e^ n (r) - a mn (r) (T - T 0 )) (28) 

is the nonisothermal increment in strain. The tensors AD^kfiv) and Aaw(r) represent the 
incremental changes in the elasticity and thermal expansion tensors due to the temperature 
increment AT. 

In a unified viscoplastic constitutive formulation [34] which is integrated by an explicit 
forward difference method, the inelastic strain increment Ae kl (r) is a function of the current 
stress (at the beginning of the increment), <r,j(r), and the current values of the internal 
viscoplastic state variables, qfi r). For example, if 

£ij ~ fij ( a rs, Qs) (29) 

then Asfj = f tJ (cr rs , q s ) At, and the inelastic strain increment is independent of the total 
strain increment Ae kl (r). This independence of the inelastic strain increment on the total 
strain increment is no longer true if an implicit integration method ( e.g . backward difference) 
or subincrementation method is used. 

The elasticity tensor D, /W (r) may be written as 

D iikl (r) = D;; tl + 6D ijt M (30) 
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where 


«W r) = 0(r) (Dl m - D%,) (31) 

and $(r) = 1 in the fiber and d(r) = 0 in the matrix, the superscripts / and m referring to 
the elasticity tensor of the fiber and matrix, respectively. The constitutive equation at any 
point r can then be written, from equation 26, as 

Aay(r) = r)) (AeJ, + Ae k ,(r) - Ac u (r)) (32) 

or 

= D? jU (Ae% + A £w (r)) - 

- (D," a Acii(rl - 6D ijU ( r) [AeJ, + Ae u (r) - Ac w (r)] } (33) 

If the quantity in braces is set equal to D™ kl Ae kl (r), that is, if 

D ijkA4i ( r ) = D ijkA c ki(r) - 6D ijU ( r) [Aeg, + Ae fc ,(r) - Ac kl (r)] (34) 

then equation 26 can be written in the form, 

Arrjj(r) = D 7 ^ kl (Aejy(r) — As£ z (r)) = D™ kl (Ae]y + Ae k i(r) — Ae* kl {r )) (35) 


From the preceding equation it is evident that the eigenstrain increment, Ae^(r), represents 
the incremental deviation from isothermal elastic behavior in the composite material when 
the elasticity tensor is taken to be a spatially constant tensor appropriate to that of the 
matrix phase. 


Newton’s law for continuing static equilibrium throughout the 

strain increment requires 

that 


dAvtjir)) 

dxj 

(36) 

Equations 35 and 36 then require that 


a{^,(^ + A £ «(r)-A4,(r))} =o 

OXj 

(37) 

or, if Ae° kl is constant, 


n m d(A£ k i(r)) d(Ae* kl (r)) 

ijkl d Xj ijkl dxj 

(38) 


4 Fourier Series Approach 

4.1 Fourier Expansions 

The application of Fourier series to the calculation of the “effective” overall constitutive 
behavior of periodic composites has been dealt with in detail by Nemat-Nasser and his col- 
leagues [25,26,27,28]. This work is used in this section to develop constitutive relationships 
for viscoplastic composite materials under small displacement conditions. 


10 


Due to the geometric periodicity of the composite we may expand Au k {r) and Ae* kl (r) in 
a Fourier series (c/., for example, Appendix 3 of Mura’s book, [24]). This gives 


±00 ±00 ±00 (2tTTI\ 27T712 21X71$ \ 

Au t (r)=£ Y. E' ^(n l ,n 2 ,n 3 )ei~ n *~^~^ n ) 


(39) 


n\ =0 n 2=0 ri3=0 


where L\, L 2 , L 3 are the dimensions of a unit cell in the x\, x 2 , x 3 directions. The coefficients 
A iik in the Fourier expansion are determined by multiplying each side of equation 39 by 

' 2nm\ 27X7H2 2ixm$ 


( Itx m 
~ 


' X 1 f X2 "b r X 3 I 

L2 Ls ‘ and integrating over the volume of the unit cell to give 


1 ^ 1 L>2 L$ , 

A u k (ni,n 2 ,n 3 ) = III Aujfc ( r ) e * 


( 2 ixn\ 2 txu2 2 txti$ 


Li Xi + ~I^ X2+ ~lT X3 


) dx 1 


dx 2 dx 3 (40) 


Xi=0 X2=0 X3=0 


where only the terms with m* = rii survive in the summations. 
Equations 39 and 40 can be written in shortened form as 


±00 


Au fc (r) = (£) e * r 

Tip— 0 

with coefficients A u k (£) determined by the inverse relation 

A u k (£) = ~r JJJ Au k (r)e~^ T dV (r) 

c v c 

where 


with 


£ = (£1 >&»&)> r = {xi,x 2 ,xz ) , V e = LxL 2 L 3 

27771; 


Si = 


L; 


(no sum on i) for i = 1,2, 3. 


The strain increment Ae* kl (r) can also be expanded in a Fourier series to give 

±00 


Ae;,(r) = 


(41) 


(42) 

(43) 

(44) 

(45) 




with coefficients A i kl determined by the inverse relation 


Ai;, (€) = ~ JjJ Ael,(T)e-+' dV(T) (46) 

C Vc 

In equations 41 and 45 the prime indicates that the term with rii = n 2 = = 0 is excluded 

from the summations, since Au*, (rii = 0, n 2 — 0, n 3 = 0) represents a rigid body displacement 
increment and A § kl ( n\ =0,n 2 = 0,n3 = 0) represents a spatially uniform strain increment. 
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4.2 Equilibrium Equation 

By substituting equation 41 into equation 19; equation 19 into the left hand side of equation 
38; and equation 45 into the right hand side of equation 38, the equilibrium relationship 
becomes 

±00 1 ✓ X 

Di‘ u EEE'»( 4 ”t (€) 00 + A*, (() hi, ) r = 

n p = 0 Z V 7 

±oo 
Tip — 0 

or 

D^hiAuk ({) = -iDT^iASi, (() (47) 

If £ = denotes the magnitude of the vector £, a unit vector £ in the direction of £ can 
be written as Q = £*/£. Equation 47 can therefore be written in the form, 

i 2 D™u (|) (|) AM«) = -iD™ a iAi' u {t) 

or 

i 2 (d^QQ) Au, «) = -«£>' J u iAti,(() (48) 

The second rank tensor, 

Mil, «) = M u (<) = D™ U Q i, (49) 

is called the Christoffel stiffness tensor (c/. [33]) and equation 48 can be written as 

i 2 M, t (<) Au k (S) = -iD^iAK, (€) (50) 

This equation can be inverted by premultiplying each side by the inverse tensor £ -2 M -1 to 
give the Fourier expansion coefficients 

(*) = -iMr k l (C) DT jr AAKs (0 T 2 (51) 

The expansion coefficients can now be substituted into the Fourier expansion of Au^r) in 
equation 41 to give 

±oo 

Au t (r) = - £ Y. T! iC 2 M- k l (0 DT ir ,iAl'r . «) e* P (52) 

Tip —Q 

This result may now be substituted into equation 19, so that the perturbation strain increment 
may be written as 

ioo 1 / v 

Ae«(r) = EEE'j ( C 2 M,;' (c )iii, + i~ 2 M a ~' (Oijh) DT jr ,Ai'„ (€) r (53) 

n p =0 Z V ' 

If we define the fourth rank tensor g k uj (£) by the relation 

m, (0 = \ (m£' (0 00 + M (0 00=) (54) 
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( 55 ) 


then the perturbation strain increment can be written in the form 

±oo 

Ae«(r) = E E E Sun (0 D? jr ,te'r, («) e if " 

Tip — 0 

and by inserting the relation for the Fourier expansion coefficients A i* a from equation 46, we 
obtain 

-I ±oo 

Ae fc/ (r) = - 

Vc n p = 0 

where the integration extends over the volume, V c = L\L^L^, of the unit periodic cell. 

From equation 18 the total strain increment is given by 

AeS(r) = A4 + 1 E E Y,' Sun «) JJJ (r') ^ dV( r') (57) 

C n P=° J Vc 

which, from the definition of D^ ra Ae* s (r') in equation 34, may be written in the final form, 

A4(r) = A4 + y E E Y'sun «) fjh ,ur - r ‘ ) {DZrACr. 00 - 

C n p = 0 y c J 

- 6D ijr . (r') [Aejr. (r') - Ac„ (r')] } dV( r') (58) 

This implicit integral equation — equation 4 in section 2.1 — must be solved to yield the total 
strain increment Ae^r) at each point r in the unit periodic cell. 

Instead of solving for Ae kl (r) from this implicit integral equation, we could use equation 34 
to eliminate Ae kl (r) from equation 57 to give an equivalent integral equation for Ae* k[ (r), viz., 

= DT ik Acu( r) - SD ljU ( r) [AeJ, - Ac ti (r)] - 

1 -Loo r c C 

- 6D ijU ( T) v EE E' sunn (0 jjj DZnrAK, (■*') dV(r') (59) 

c n P=° Vc 

The incremental constitutive relation at any point r is given in equation 35, and this 
relation can be used to update the stress state at any point r in the unit cell once equation 59 
is solved for Ae^(r). Alternatively, equation 58 can be solved for Ae kl (r) and inserted into 
equations 34 and 35. The overall “effective” constitutive relation for the composite material 
can be obtained by averaging equation 35 over the unit periodic cell. This gives 

(A <r„> = (A<4 + A e u - Ael,)) 

or 

(A <t„) = D"‘ t A4, + <Ae«) - (A 4> (60) 

If we equate the volume averaged stress increment (A er.jj) and the overall bulk response stress 
increment A of-, i.e., if (A o { j) = Act? , and we note from equation 24 that the volume averaged 
perturbation strain increment is zero, i.e. (Asm) = 0, then the overall “effective” constitutive 
relationship is 

A<4 = D? jU Ae», - D^u (A 4) (61) 


K) III DTjrAK, 00 e it (r - r '> dV(r') (56) 

V c 
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or 


A<7“ = D% t ,Ae° u - i JJJ {DT SU Ac*, (r) - SD«« (r) [Ae£ (r) - Ac*, (r)] } dV (r) (62) 

° Vc 

which is the result presented in equation 2 of section 2.2. The procedure for integrating the 
overall “effective” constitutive relation then proceeds as follows. 

4.3 Fourier Integration Algorithm 

1. From a knowledge of the stress state throughout the unit periodic cell at 
the current time, t, calculate the inelastic strain increment Aej^ (<r rs , r) from an 
appropriate unified viscoplastic constitutive relation. The viscoplastic constitutive 
relation will vary according as r is in the fiber or matrix phase, respectively. 

2. Compute the eigenstrain Ae^(r) throughout the unit periodic cell from the 
implicit integral equation 59 or from equations 34 and 58. 

3. Compute the stress increment throughout the unit periodic cell from equa- 
tion 35 and update the stress, strain and viscoplastic state variables according to 
the relations 

® ij (r, t T At) o ij (r, t) T Aay(r), 
efj ( r > t + At) = efj (r, t) + AeJ(r), 

Qi (r, t + At) = qi (r, t) + A q t (r) . 

4. Calculate the overall “effective” stress and strain increment for the compos- 
ite from equation 61 and update the overall “effective” stress and strain from the 
relations 

a% (t + At) = 4 (t) + A<r» , 

£ ij (t + At) = (t) + Ae° . 

5. Repeat the preceding calculations for each incremental load step. 

4.4 Implicit Integration Algorithm 

The preceding algorithm makes use of the fact that the inelastic strain increment Ae^(r) is 
independent of the total strain increment Ae^(r) if an explicit forward difference method — 
such as Euler or Heun forward difference — is used to integrate the unified viscoplastic relations 
for the fiber and matrix phases. If an implicit method — such as backward difference or sub- 
incrementation — is used, the inelastic strain increment depends on the total strain increment. 
In this case the total strain increment must be obtained by iterating equation 58 in the form, 

Ae[,(r) = A4 + y Y. E Y!sm (0 /// ^ ^ { D ?jrA c r, (r', AeJ, (r')) - 

c n p =0 y- 

- 6D IJr , (r') [a 4 (r') - A c,., (r\ AeJ, (r'))J } dV(r') (63) 
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The first iterative guess can be taken as Ae^(r) = Ae°;, and the right hand side evaluated 
to give an improved guess for Ae kl (r). This process is then continued with 



f E £ E'itey (0 /// <rr,> {D™,Ac„ (r', {A e £ (r')}J - 

c n p = 0 v/ 

- (O [{A^ (r')} A - Ac rs (V, {AeJ, (r')} J] } dV{y') (64) 


until the A th and (A + l) th iterates of Ae kl (r) converge. 

Equation 59 is not so convenient for iteration as equation 58 when Ae^(r) depends on 
Ae kl (r). It is always necessary to know the total strain increment Ae kl (r) in order to calculate 
the inelastic strain increment Ae^ s (r', Ae^ (r')j. But equation 34, viz., 


D m A? 

L Ujk i^ £ 


«M = D ?jkA c “ (r, AsJ/rj) - 6D ijkl ( r) [Ae£(r) - Ac« (r, A £ J,(r)) 


(65) 


is an implicit equation for Ae kl (r) when the iterated quantity, Aej^(r), is given. Equation 63 
is therefore the appropriate equation to iterate when the inelastic strain increment depends 
on the total strain increment. 

The procedure for solving the implicit integral equations in 58, 59 and 63 is described in 
section 8. 


5 Green’s Function Approach 

5.1 Green’s Solution of Navier’s Equilibrium Equation 

The equation of continuing static equilibrium for the composite material throughout an ap- 
plied strain increment is given by 


a(A ^ ,(r)) + A/j(r) = 0 (66) 

where A/, (r) is the incremental body force per unit volume of the composite material. From 
equations 35 and 66 we obtain 



which is equivalent to equation 37 in the absence of the incremental body force A/)(r). From 
this equation it is clear that the divergence of the stress variation produced by Ae^(r) may be 
formally regarded as a fictitious body force increment, analogous to A/,(r), which is applied 
to the homogeneous matrix material with elasticity tensor D™ kl . The theory of elasticity for 
homogeneous materials is generally concerned with the solution of the homogeneous differ- 
ential equation 67- Navier’s equation — when the right hand side is zero. When body forces 
are present the standard method of solution is to obtain the displacement solution at r due 
to a unit body force applied at r'. This solution is given by the Green’s function (r — r') 
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which gives the displacement in the i th direction at r due to a unit point force applied in 
the j Ul direction at r'. For a distributed incremental body force A/, (r') the displacement 
increment at r is obtained by summing the results for the distribution in the form 


Altj(r) = fjfc,, (r - r')A/, (r') dV( r') (68) 

The integration extends over the whole volume, V . of the composite material which may be 
regarded as being of infinite extent. 

When = 0 we know that the displacement solution is Auf(r) = Au°(r), cor- 

responding to an applied uniform strain increment Ae° on the infinite boundary of the 
homogeneous matrix. For an effective distributed body force increment given by the right 
hand side of equation 67, with A fj (r') = 0, the solution for the total displacement increment 
Auf(r) can be written as 

AuJ(r) = Au“(r) - (r - r') A AeL.(r')) dV(r') (69) 

V 1 

This corresponds to equations 17 and 39, the volume integral corresponding to the perturbed 
displacement increment Auj(r) in 17. 

For a material which is homogeneous with elasticity tensor D™ kl the Green’s function 
satisfies the differential relation (c/. Appendix A, equation A. 11), 

nm 9 2 Gi bn(r-r') * ^ n 

dx.Bx, + (r - r ) = 0 (70) 

where 6 im is the Kronecker delta tensor given by S im = 1 if i = m and 6 im = 0 if i ^ m, and 

8 (r — r') is the three dimensional Dirac delta function defined by the relation 

6 (r - r') = 8 (aq - x\) 6 (x 2 - x' 2 )6(x 3 - x 3 ) (71) 


By applying the Fourier integral techniques in Appendix A, the Green’s tensor is shown to 
have the Fourier integral form, 


(r - r>) = /// 


d 3 K 

(2tt) 3 


M tj 


-l 


K 2 




iK.(r-r') 


(72) 


in which the inverse Christoffel stiffness tensor (c/. [33]) M tJ 1 (£) is defined by 

M ,- 1 «) = < 73 > 

with ( p = K p j sj K m K m — K p /K being a unit vector in the direction of the Fourier wave 

vector K, and K = JK m K m denoting the magnitude of the wave vector K. 

Making use of the relation 

G lt (r-r')A ( DSmnA<m{r '^ = A ( Gji (r - (r')) - 
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we may write equation 69 in the form 

Auf(r) = Au°(r) - JJJ JL (a* ( r - (r')) dV( r') + 

+ Iff 9G,k £r~ D? lmn iu' mn (r') dK(r') (75) 

V z 

The first volume integral can be transformed into a surface integral via Gauss’ divergence 
theorem, viz., 

JJJ { Gik (r “ r ^ D kimn^*mn W)) dV (r') = JJ n, (r') G ik (r - r')D% mn Ae* mn (r') dS( r') 

V 1 s 

( 76 ) 

The surface integral extends over the entire outer surface of the “infinite” matrix material. 
Since this is assumed to be at an infinite distance, all the integration points r' in the surface 
integral are at an infinite distance from the field point r and G ik (r — r') = 0. Thus, for an 
infinite body the first volume integral in equation 75 vanishes. This would not be the case 
for a finite body in which the field point r is close to the surface integration point r', and 
the volume (or surface) integral would need to be retained for these situations. In this case 
other surface integrals would arise (c/. Appendix D, equation D.27) due to the application of 
boundary incremental displacements or surface tractions on the surface of the material. 
From the properties of the Green’s function, 

dGg fc (r - r') _ dGjk (r - r r ) 
dx\ dxi 

which follows since G ik is a function of 

r - r' = (xi — x[,x 2 - x' 2 ,x 3 - x' 3 ) 

Equation 75 may then be written alternatively as 

A uj (r) = A«“(r) - Jj/ ^^ ~ ^- DZ mn Ae‘ mn (r') dV( r’) (79) 

v 1 

But AsJ(r) = %(d (A uf (r)) /dxj + d (A uj (r)) /dx l ), so that by differentiating equation 79 
with respect to Xi and Xj and taking half the sum, we obtain 

AeJ(r) = A + JJJ U m (r - r') D% mn Ae* mn (r') dV (r') (80) 

v 

which, by means of equation 34, may be written as 

A 4( r ) = A 4 + JJJ U m (r - r') {D% rs Acr S (r ; ) - 

- 6Du „ (r') [AeJ, (r') - Ac„ (r')] } dV(r') (81) 


(77) 

(78) 
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An equivalent integral equation, involving the eigenstrain increment Ae^r), can also be 
obtained by using equation 34 to eliminate Ae^(r) from equation 80, which gives 


DljkA4l M = A7w^c«(r) - SD m ( r) [AeJ, - Ac u (r)] - 

- SD m ( r) JJf Ukimn (r - r') DZ„A<, (r') dV( r') 


(82) 


In the preceding equations the operator, 


Uijki (r - r') 


1 l d 2 Gifc (r - r') (r - r') \ 

2 \ dxjdxi dxidxi ) 


(83) 


gives the ij component of the strain increment at point r due to an applied stress increment 
component kl at point r' in an infinite homogeneous medium with elasticity tensor D”J kl and 
Green’s function given by equation 72. 


5.2 Equivalence of Perturbed Strain Increment 

From equations 18, 56 and 80 we see that the perturbed strain increment, Ae w (r) = Ae kl (r) — 
Ae k/ , is given by the equivalent relations, 


"I ioo c r c 

Ae«(r) = 77 E £ Y^dkimn (0 JJf D™ nrs Ae*rs (r') e dV( r') (84) 

Vc n P=° Vc 

or 

AE “ (r) = JfJ Uu (r - r '> D ™r* A <‘ M dv (r') (85) 

V 

The volume integral in the Fourier series representation extends over the volume, V c , of 
the unit periodic cell and the summation extends over the integers n p = 0, ±1, ±2, . . . , etc., 
where p = 1,2,3. In the Green’s function approach the volume integral extends over the 
entire infinite medium, i.e., over all the periodic cells comprising the material. It is shown in 
Appendices B and C that the Fourier summation expression in equation 84 can be converted 
into the Green’s function expression in equation 85 by means of the Poisson sum formula. 

From equation 34 it is evident that if the elastic properties of the fiber are the same as 
that of the matrix, then 8D,jki{ r) = i9(r) (Df jkl — = 0, in which case 

A 4/( r ) = Ac w (r) (86) 

is known explicitly without having to solve the integral equation. From equations 58 and 81 it 
can also be observed that Aejr(r) is known explicitly when 8Dijki( r) = 0. The explicit relation 
in equation 86 holds only when an explicit forward difference method is used to integrate the 
viscoplastic constitutive relations. For implicit integration methods in which the inelastic 
strain increment Ae^(r) depends on the total strain increment Ae^(r), equations 58 and 
81 show that even when <*> 7\,«(r) = 0, the equation to determine As kl (r) is still an implicit 
integral equation. 
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6 Self-Consistent Method 

6.1 Outline of Self-Consistent Method 


In this section we establish a self-consistent relationship between the overall “effective” stress 
increment, A<r° , and the applied strain increment, Ae° , for a matrix material which has 
cylindrical fibers embedded in it in a periodic fashion. 

From equations 34 and 61, this relationship can be written as 

= DT tU A4 - ^jjj {D™ kl Ac tl (r) - SD m (r) [Ae£ (r) - Ac„ (r)] } dV( r) (87) 

C JJJ 

where the total strain increment is determined from equation 81 in the form, 

A4(r) = Ae° kl + JJJ U klmn (r - r') {£>^ nrs Ac rs (r') - 

V 

- 6D mnrs (r') [AeJ s (r') - Ac rs (r')] } dV ( r ') (88) 

These equations can be solved in an approximate fashion by means of a self-consistent 
method in the following manner. 

First, assume that the unit periodic cell consisting of a cylindrical fiber embedded in a unit 
matrix cell, Fig. 5, is replaced by a cylindrical fiber (of radius = a) embedded in a cylindrical 
matrix (of “effective” radius = 6) as depicted in Fig. 6. The other unit cells outside the given 
unit cell — i.e., the rest of the composite — are then smeared out into a uniform matrix material 
whose overall “effective” constitutive properties are the volume average of the constitutive 
properties of the constrained unit periodic cell. The “effective” constitutive properties will 
be transversely isotropic if the fibers are arranged in hexagonal arrays or tetragonal if they 
are arranged in square arrays. 

Second, assume that the total strain increment, Ae kl (r), and the strain increment repre- 
senting the deviation from isothermal elastic behavior, Ac^(r), are spatially constant in the 
fiber and matrix phases of the unit cell. These constant values (different in the fiber and 
matrix of the unit cell) are taken to be the volume averages over the respective constituent 
volumes of the fiber and matrix phases of the unit cell. 

The composite now consists of three constituent phases, viz., the fiber, matrix, and 
smeared out average phases. If the elasticity tensors of these phases are denoted by Df jkl , 
D™ kl and D ijk i, respectively, then the elasticity tensor at any point r in the composite can be 


written as 

D k lrs{^) D k lrs T ^Dklrsi, r) 

(89) 

where 

fiDklrsi r) — 8D) hrs — D k i ra — Dklrs 

(90) 

if r is in the fiber; 

bD k i TS {r) = 8D' klrs — D 7 klrs — D k i rs 

(91) 

if r is in the matrix; and 

SDkirsi r) = 0 

(92) 
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if r is in the surrounding smeared out “effective” material. 

The fiber and matrix constituent phases now represent fictitious body forces in the infinite 
“effective” medium with elasticity tensor Dijki , and the total strain increment is obtained from 
the solution of the integral equation, 


in which 


A4(r) = Ae% + JJJ Uiju (r - r')D tlr ,Ae‘, (r') dV( r') 

V 


D klrsA£ TS (*" ) — ^klraAc rs (f ) $Dklrs (*" ) 


(r') - Ac rs (r') 


(93) 

(94) 


6.2 Strain Increments in Three Phases 

We now make the approximation that the strain increments in the three phases are spatially 
constant and equal to their respective volume averages, so that if r' is in the fiber Aefj (r') 
and Acij (r') are replaced by 

Ael(f) = yjjjAel(Y')dV(T') (95) 

f Vf 

and 

Ac H (/) = Y JJJ Acy (r') dV(r') (96) 

5 v f 

so that, from equation 27, 

Acij(f) = Asfj(f) + AT (97) 

with 

a i/ AT = q{jAT + (T — To) Aa{j — 

- (D{ jki y' AD' lmn (eL </) - <&,(/) - <*L ( T - To)) (98) 

If r' is in the matrix the relations are replaced by 

Ae? (m) = 77 “ JJJ Ael (r') dV( r') (99) 

m 

v tn 

and 

Ac tJ (m) = Jj [Ac,, (r') dV( r') (100) 

m 

v m 

where 

Acij(m) = A e?(m) + a)J"AT (101) 

and 

ay" A T = n"jAT + (7 1 - T ( >) Ao'" - 

- (D^,y [ ADH m „ (4„(m) - eL(m) ~ <n W ~ To)) (102) 
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If r' is in the smeared out “effective” or homogenized medium the corresponding results are 



A 4 = A 4 = v III <r,) dV(r,) 

s v/. 

(103) 

and 

Ecy = V III A ° ij ^ dV ^ 

s v a 



(104) 

where 




Ar Cij = Ai^ + a'AT 

(105) 

and 




a* AT = ciijAT + (T — T 0 ) Aa^ — 

— ^ijkl ADkhnn mn ~ £ mn ~ a mn (T ~ ^o)) (10®) 

The volumes V), V m and V s refer to the volumes of the fiber, matrix and smeared out medium, 
respectively. If V c is the volume of the unit cell and V denotes the total volume of the entire 
composite material, then 


V c = Vf + V m and v = V c + V s = Vf + V m + V s (107) 

6.3 Applied, Homogenized and Volume Averaged Increments 

At this point it is important to emphasize the following distinctions. First, the strain incre- 
ment applied to the composite is denoted by A which causes an incremental stress response 
A < 7 ° . To obtain the overall “effective” constitutive equation these are equated to the cor- 
responding volume averaged quantities, ^A and (Aer^). In the “effective” homogenized 
medium all quantities are denoted with overbars. 

At any point r the appropriate constitutive relation is 

Acr^r) = D m (r) (Ae^r) - Ac fc/ (r)) (108) 

If we volume average this relation over the unit cell we obtain 

(A c?ij) = (^DijkiAel^j — ( DijkiAcki ) (109) 

In the homogenized phase the constitutive relation can be written as 

Act ij DijkiAe kl D^j^iAc^i ( 110 ) 

T 

Since the strain increment Ae kl in the homogenized phase must correspond to the applied 
strain increment Ae kl ~ as in equation 103 — and the homogenized stress increment Acr^ must 
correspond to the overall bulk stress increment Act?-, we write the constitutive relation for 
the homogenized phase as 

Aofj = DijkiAch - Dij^Acki (HI) 


21 



6.4 Requirement for Self-Consistency 

For self-consistency we require that the volume average of the microscopic constitutive relation 
in equation 108 over the unit cell, viz. equation 109, should correspond to the constitutive 
relation for the overall “effective” homogenized medium in equation 111. That is, 

ACjj = (Dijki&Cki) — DijkiAEfci DijkiACkl (112) 

for self-consistency. Under the approximation that the strain increments Aeh and Ac*.* are 
spatially constant in the constituent phases, we obtain 


A<r“. = A /// D, jtl (r) (A4(r) - Ac M (r)) dV( r) = D m AeJ, - D m Ac« (113) 

v c 


or 

A<T?j = y D l 3 ki (teL(f) - Ac,w(/)) + ( a 4(™) - &c k i{mj) 

= DjjkiAs DijkiAcu ( 114 ) 

At this point the elasticity tensor Dijki and the deviation strain increment Acki in the homoge- 
nized medium are unknown quantities. In the next section we will solve the integral equations 
for the total strain increments in the fiber and matrix phases, Aeh(f) and As^(m), and we 
will find that these values depend on the quantities D^i, As jy and Acf-i in the surrounding 
homogenized medium. Then, by equating the coefficients of A e Q kl on both sides of equa- 
tion 114 we obtain a relationship for the unknown elasticity tensor Dijki of the “effective” 
homogenized medium. The value of the unknown deviation strain increment A Cki in the ho- 
mogenized medium can then be obtained by equating the terms independent of Aejy on the 
left hand side of equation 114 to the corresponding term DijkiAcki on the right hand side. 

We now obtain the total strain increments Ae^(/) and Ae^(m) in the two phases of the 
unit cell. First, consider the total strain increment in the fiber phase. 


6.5 Total Strain Increment in Fiber Phase 

Equation 93 can be volume averaged over the fiber phase to give 

A f[f tel(T)dV(r) = A 4 + A JfJdV(r)ffJ, U ijU (r - t')D kl „U' r , (V) dV( r') (115) 

/, v> f Vf v 

where the field points r are in the fiber volume, Vf, and the integration points r' are in all 
three volume phases {V — Vj + V 1U + U s ). Equation 115 can be written as 


Aefj (/) — A £®j + — 


JJJ dV( r) JJJu m (r - r') D u „ As’. (/) dV(r') 

. Vj Vf 

+ f f j U,j kl (r - r') Dkirs Ae* rs ( m ) dV( r') + JJJ U ijkl (r - r') D k i rs Ac rs dV{ r') 

' i';„ i/ s 


(116) 


22 



(117) 


in which Dkin, Ac* s (r') has been replaced by 

D Ur AK, (/) = D Ur ,Ac r ,(f) - 6D' Ur , [A e T „(f) - A .*.(/)] 

and 

D k i rs Ae* a (m) = D klrs Ac rs (m ) - 6D% rs [Ae£(m) - Ac^m)] (118) 

in the respective fiber and matrix phases, and by 

D k i r3 Ae* s (s) = D k i rs Acr S (119) 

in the smeared out “effective” medium where 8D k i rs (r') = 0. 

In the first integral in equation 116 the field point r lies in the volume Vf, and since 

JIJ U ijU (r - r') dV(r') D klrs = S ijrs (120) 

v 

is Eshelby’s tensor (c/. Appendix E, equation E.l and [35]), which is a constant tensor 
independent of r when the field point r lies within the cylindrical volume V included in an 
infinite medium with elasticity tensor D k i rs , we may write the first integral as 

JJJu tjkl (r - r') dV(r') D klrs A£* TS (f) = S ljrs Ae* rs (f) (121) 

Vf 

The second volume integral extends over the volume V m — V c — Vf of the matrix phase. Thus, 
for the second integral, 


JJJ U m (r - r) dV( r') D klrs AE* rs {m) 

Vm 

= fffu m (r - r') dV(r') D klrs Ae* rs (m) - Jjju m (r - r') dV(r') D klrs Ae* rs (m) 

Vc Vf 

= S ijrs A£* rs (m) - Sijrs A£* rs (m) = 0 (122) 

since the field point r lies in the cylindrical volume Vf and therefore within the cylindrical 
volume V c . 

We now have to deal with the last integral in equation 116. This integral can be written 
as 


JJJu.jU (r - r') dV( r') £> W) , s Ac„, 

Vs 

[[[- ( 9 d^tr-r') d 8G jk (r - r') 

JJJ 2 \dxj dxi dxi dxi 


dV (r ) D k i rs Ac rs 


fffl ( d dGik{ r-rQ d dG jk (r - E) 

JJJ 2 l dx\ dxi dx'i dxi 

\ r s J 


dV{ r') D k i rs Ac rs 
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which can be transformed via Gauss’ divergence theorem into two surface integrals: one over 
the outer “infinite” surface S of the smeared out “effective” medium; and the other over the 
inner surface of the “effective” medium, i.e., over the surface S c of the unit periodic cell. The 
volume integral then takes the form, 


JJJ Uijki (r - r') dV(r') D k i rs Ac rs 
v. 


= // 1 (»; M - *1^ + n* <r') dG * (r - r ' ) ' 


+ H\(w 

Sc V 


dxi 

dG ik (r - r') 


dxi 


dxi j 


dS(r') D klrs Acr S + 


dS(r ) D k i rs Ac rs 


where n* (r') is a unit normal vector at point r ; on the surfaces pointing away from the volume 
V s . Now since the field point r lies in the fiber and the integration points r' on the surface S 
are infinitely removed, we have dG ik (r — r')/dx{ — > 0 on the outer surface S of the composite, 
and the first surface integral can be neglected. If we write n, (r') = — n* (r'), then n, (r') is a 
unit normal vector pointing away from the volume V c on the surface S c of the unit cell, and 
we have, via Gauss’ divergence theorem and equations 77 and 83, 


1 ( d dG ik (r — r') , d dG jk ( r - r') 


JJJ U ijkl (r - r') dV(r') D klrs Ac r 

Vs 

ff 1 ( , n dG ik (r - r') , 

= -JJ 2p (r) ^ +n ' (r) 

Sc V 

* 'Jill 

-k 

Vc 

= - JJJ U ijk i (r - r') dV( r') D klrs ~Ac r 

Vc 


.. dG jk (r-r ')' 


dx t 


dS(r ) D k i rs Ac rs 


dx', 


dxi 


- + 


dx'c 


dx t 


dV( r') D klrs Ac r 


1 ( d 2 G ik (r - r') d 2 G jk (r - r') 


dxjdxi 


+ 


dxjdxi 


dV (r ) D k i rs Ac TS 


Since the field point r lies in the cylindrical volume V c , the preceding equation takes the form, 


or 


JJJu ijkl (r - r') dV( r') D klrs Ac rs = - JJJ U m (r - r') dV( r') D klrs Ac r 

Vs vj 

J ' U{jki ^ ) dV(r ) DklrsAcrs Sij r $Ac rs 


where 5, jr , s is Eshelby’s tensor for a cylinder with elasticity tensor Dij k i- 
From equations 116, 121, 122 and 123 we obtain 

AeJj(f) = Ae» + jt/// dV(r) 


/ f ./ J 

V> 


( 123 ) 
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and since the quantity within the braces is a constant tensor, 

Ae J(/) = Ae», + S, Jr , (Ae;,(/) - Acts) 

Now 

D klrs Ae* ra (/) = D klrs Ac rs (f ) - 6D f klrs [A £ J s (f) - A c rs (f)} 

so that 

Ae;, (/) = A Cr.U) - [AeL(/) - A<w,(/)] 

and on substitution into equation 124 we obtain, 


(124) 


A £$(/) = Ae% + S ijrs yAcra(f) - Ac rs - [A £„(/) " Ac mn (/) 

Given that 

lijkl = 2 i^ik^jl 4 " dj k dil) 

denotes the fourth rank identity tensor for symmetric second rank tensors, the preceding 
equation can be written as 

\jijmn 4“ SijrsD rspq SD^ qvm A £ mn (f) — A£- ~h 

4“ Sijrs ^Ac rs (y*) Ac r<s ^ H" Sij rs D rS pq <5Z)^ mn ,Ac y7m ( i /) 
which, by premultiplying each side with the inverse of the tensor in square brackets, gives 


A eUf) = 


^ij mn ~t“ $ijr. rspq ^ ^ pqmn 


{a £ m n Smnrs (A c rs (f) - A c rs ) + 

4- S mnk iD klpq SD^Acr.lf)} (125) 

The phase volume averaged stress increment in the fiber is then given by the relation 

A <r«(/) = D{ jU (A 4,(1) - A c„(/)) (126) 

6.6 Total Strain Increment in Matrix Phase 

Now consider a field point r in the matrix phase. From equation 93 we may write 

AeJ(r) = Ae% + JJJ U m (r - r') dV(r')D klrs Ae* s (f) + 

v r 

+ ///C„«(r- r') d,V{r') D klr sAe* ra (m) + JJJ U tjkl (r - r') dV{r') D klrs Ac rs (127) 

I'm V, 

Since V,„ — V c — Vf, the second integral can be written as 

// / u, jtl (r - r') rfV(r') D Ur , A £ ;,(m) 

' V,n 

= fJJUijkd r-r') dV{v')D klra Ae* ra {m) - JJJu ijkl (r-r') dV(r')D klrs Ae* TS {m) (128) 
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and from equation 123 the last integral in equation 127 may be written as 


Uf Uijki (r - r') dV (r') D k i rs Ac rs = -JJJ U ijM (r - r') dV (r') D k i rs Ac ra (129) 

' v,' v c 

so that equation 127 is transformed into 

AeJ(r) = Ae% + JJJu ijtl (r - r') dV( rj D Ur , (A e'„(f) - A e’Jm)) + 

V f 

+ JJJUijki (r - r') dV{v’)D klrs (A e* rs (m) - A c ra ) (130) 

y c 

Averaging equation 130 over the matrix phase gives 

A ejjm) = Ae% + 2-JJJdV(r)UjJ U ijt , (r - r') dV( rj D Ur , (A e‘„(f) - A e' r ,(m)) + 

m l Vf 

JJJ U ijk i (r - r') dV(r') D k[rs (A ^(m) - Ac ra ) 


+ 


or, since V m = V c - V f) 


(131) 


AeJ(m) = Ae% + ~ 

V m. 


JJJdV( r) JJJ Uijki (r - r') dV( r') D k i rs (A e r %(/) - A£ r * s (m)) 

Vc V/ 

+ JJJ dV( r) JJJ U^ki (r - r') dF(r') D fc/rs (A e* s (m) - Ac rs ) - 

Vc Vc 

- III dV ^ III Uljkl ^ dV ( r ')Dur S (A £* ra (f) - A£* s (m)) - 

'/ V/ 

-JJJ dV ^III Uijkl ~~ dV ( r '} Dkirs (a £* rs (m) - Ac r . s .) 


+ 


(132) 


Now consider the first integral in equation 132. We may interchange the order of the 
volume integrals so that 


JJJ dV(r) JJJ U im (r - r') dV(vj = JJJdV( rj JJJ U ijU (r - r') dV( r) (133) 

V,' Vf Vf V c 

Now r and r' are dummy integration variables, so that on the right hand side of equation 133 
the variables may be replaced with the integration variables x and y, viz., 

JJJ dV(r) JJJ U !jU (r - r') dV(rj = JJJdV( y) JJJu ijU (x - y) dV(x) (134) 

'vj ’ Vf Vf Vc 
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But, from equation F.5 of Appendix F, 


Uijki (x - y) = Uijk, (y - x) 

so that 

JJJdV(r)JJJlI m (r - r') dV( r') = JJJdV( y) fjju m (y - x) dV(x) (135) 

' Vc. Vf Vf V c 

and the dummy variables y and x may be replaced with the variables r and r' to give 

JJJdV (r) jjju iiu (r - r') dV( r') = JfJdV( r) JJJ U m (r - r') dV( r') (136) 

Vc Vf Vf Vc 

This relationship is discussed in Mura’s book ([24], page 336) where it appears under the 
heading of the Tanaka-Mori theorem. 

From equation 136, the first integral in equation 132 is integrated over the field points r 
within the cylindrical volume Vf. Since these field points lie within the cylindrical integration 
volume V c , the first integral in equation 132 may be written as 

y- JJJ dV(r) JJJ U m (r - r') d,V( r') D tlr , (A e'„U) ~ A C(™)) 

m Vf Vc 

= Y JJJ S t , „ (A e;,(f) - ACM) dV(r) 

m Vf 

= (A C(/) - ACM) (137) 

* m 

In the second integral in equation 132 the field points r lie within the cylindrical volume V c 
and so the second integral may be written as 

t^/// dV( r) JJJ U m (r - r') dV(r') D klrs (A e* rs (m) - 

m Vc Vc 

= JJJSijrs (A s* rs (m) - ^rs) dV( r) 

m V c 

= £s iJr , (A e* ra (m) - Ac-rs) (138) 

* m 

In the third integral the field points r lie within the cylindrical volume Vf and so 

YJJJ dV(r)JJJu ljkl (r - r') dV( r') D Ur , (Ae*,(/) - A CM) 

"* < f Vf 

= Y .(/) - A <.w) dV M 

m Vf 

= Vj-S ijr . (Ae;,,(/) - ACM) (139) 
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Finally, in the fourth volume integral, the field points r lie in the cylindrical volume Vf. 
Since this lies within the cylindrical integration volume V c , we have 

V~ /// rfV(r ) III Uijkl {r ~ r,) dV ( r ')Vkirs (A < s (m) - A c rs ) 

"* Vf Vc 

= ~ III s ljr , (A e,*,(m) - S3,,) dV( r) 

= hls y „ (Ae r *,(m) - A3,,) (140) 

We thus obtain from equations 132, and 137 to 140, 

AeJ(m) = Ae% + ^5 i?rs (Ae* rs (f) - Ae* rs (m )) + 

I'm 

+ ^ (A< s (m) - Ac r ,) - 
~ ^.jr S (A £* rs {f) - A £* rs (m)) - 

* m 

Vf / \ 

T7~^ijrs 

I'm 

= As)) + (Ae;,(m) - S3,,) 

or 

AeJ(m) = Ae° + Si in , (A e* rs (m) - A c rs ) (141) 

This relation for the total strain increment in the matrix phase is similar to that for the fiber 
phase given in equation 124. By following the steps leading to equation 125 the expression 
for A efjim) can be put in the form 

A £ij( m ) — [-fijmn + S i]rsD rspq 6D pqmn j { Ae mn + S rnnrs ^A c rs (m) — Ac rs ) + 

+ S„„ m D^SV^ .Ac„(m)} (142) 

The phase volume average stress increment in the matrix is then given by the relation 

A <Tij(m) = D™ kl (A eh(m) - A c«(m)) (143) 

6.7 Overall “Effective” Constitutive Relation 

As stated in section 6.4, for self-consistency we require that the volume average of the con- 
strained micromechanical constitutive relation over the unit periodic cell should correspond 
to that for the “effective” homogenized medium. From equation 114 we require that 

Acr” = y^D f ijk i (Ash(f) - A c«(/)) + ~D"l kl (A e£(m) - Ac w (m)) 

= DijkiAeh — DijiciAcid ( 144 ) 
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where the total strain increments in the fiber and matrix phases are given by equations 125 


and 142 as 



^ 4(f) = 

lijmn + Sij r3 D rspq 6D pqmn J | A6 mn + S mnrs ^A C rs (f) A C rs J + 

+ SmnklD klpq ^Hpqrs^ C rs(f)} 

(145) 

and 



AeJ(m) = 

lijmn Sij rs D rS pq ^Hpqrnri\ ^mn ^ mnrs ^ACy^TH.) ACrs'j “1“ 

+ S mnkl D kl l q 6D™ qrs Ac rs (m)} 

(146) 

with the deviation strain increments defined in equations 97, 101, 105 as 



Ac ij (f) = Ae? j (f) + a *jAT 

(147) 


Acij(m) = Ae^(m) + a*™ AT 

(148) 

and 

A Cij = As** + a*jAT 

(149) 


By inserting equations 90, 91, 145, 146 147, 148 and 149 into equation 144 the relationship 
for self-consistency requires that 


ij {AijklSklrsACrs CijrsA£ rs {TYl) j” 

- {AijkiSki rs a? s — Bi jrs a*{ - C ijrs a*™} AT 


in which 


-FT P 


Dijki&.e k i Dijkiotfri&.T 


-A-ijkl — 


y Bjj rs Irskl 4" Brspq Dpq ?nn (^D mn kl 

-FT -1 


D{ r 

Vrn 

V c 


+ — d™ \i 

~ - - ^ijrs [ 


rskl 


i q n A ( n m 

' "~ >rs Pq- LJ pqinn y^mnkl 


Hmnkl^ | 


(150) 


(151) 


and 


Bijkl 


= Yldi. 

Y ijrs 


Irskl + UrspqDpqmn 


( D 


mngh 


-1 


X 


Ughkl “I” ^ ghabU a bcd cdkl ^cd/c/^J y I^ijkl 


Qjki — 


V r 


y ^ ijrs 


*rskl 


Q D 1 

u r spq - Ly pq ini i 


i D < 


mngh 


— D 


mngh 

v m 


r 


X Sghkl + SghabD abcd {D™ dkl - D cdk i ) - ~^D™ kl 


(152) 


(153) 
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These results for B tjk i and C^ki can be simplified somewhat. We write Bijki as 


B,jkl yr yijkl 


where 

and 

We then obtain 


yijkl Dijpq ( Ipqrs T X pgr a ) (Sr ski T X rs kl ) Bfj 


ijkl 


Xpqrs “ S p qm nD mngh SDg hrs 


Uijkl T D ijki D ijpq ( Ipqrs T X pqrs ) {S rskl T ^rskl) 

or i 

(^ijkl) ( yklrnn T ^klmii) ~ ( hjkl T X ijkl ) {Sklmn ~t~ ^Zmn) 

This result simplifies to 

(j-^ijkl) yklmn "f* lijmn ~ (hjki + Xijki) 1 ( Ski mn T Xklmn) 

which can be premultiplied by (I pqi j + X pqtJ ) to give 


(Ipqij + Xpqij ) yklmn — <S) 


pqmn I pqmn 


from which 


yijkl B ij p q (Ipqrs T X pqrs ) (S rs kl Iraki) 

From this result we find that B l]k i and Cijki can be written in the simplified forms 


Bijki — y Bijpq \lpqgh + S p q rs B rsmn (D^ngh D mng hj (Sghkl Ighkl ) 


and 


Vr 


B ijkl Y Djjpq \lpqgh + S pqrH D rsinn {^B mng h B mng hj (Sghkl Ighkl ) 


Vr 


(154) 

(155) 

(156) 

(157) 

(158) 

(159) 

(160) 
(161) 

(162) 

(163) 


Equating the coefficients of Ae kl in equation 150 for self-consistency then requires that 

Dtjki = A ijki (164) 

which, from equation 151, produces the implicit relation 


Bijki 


Y Bijrs [irskl T S r!l pqBpg jnn (jD mn kl B mn kl ^ j T 


V T 


T^-l 


+ i 1 '*' + - o—x 


-1 


(165) 


The value of homogenized “effective” elasticity tensor B { j k i may be obtained from this im- 
plicit relationship by iteration. Naturally, when the self-consistent method is embedded in a 
nonlinear finite element program, this iteration would be done outside of the code and the 
explicit values of D,jki would be used in the program. 


30 



For self-consistency we also require from equation 150 that 


AijpqSpqki&eu ~ Bijuteuif) - C ijk Ae[,{m) - D ijkl Ae kl (166) 


and 

■AijpqSpqltfOifci B>j k'l Cij k iAot k i Dij k [OC k i (167) 

which, by setting A i]Pq = D, Jpq , reduce to 


and 


ASjj — Dijp q {Spq rs - Ipqrs ) i^B rs kl^^kl(f) ”b 

Dijpq (Spq rs ~ Ipqrs )] (^B rs hiOtfci “h Crskl&kl ) 


— - ★ 

a u = 


(168) 

(169) 


The overall “effective” constitutive relation for the homogenized composite in equation 150 
is now easily computed. 

If a forward difference algorithm is used to evaluate the viscoplastic strain increments, 
the only implicit equation which occurs in the formulation is that for the elasticity tensor 
of the homogenized medium given in equation 165. It is, perhaps, ironic that in deriving 
the highly nonlinear viscoplastic constitutive relationship for the homogenized medium, the 
only iterative procedure required is that for the elasticity tensor. This implicit elasticity 
relationship also occurs in the subvolume method due to the occurrence of the tensor SD ljk i 
in the volume integration. The implicit nature of Dij k i is due to the fact that the homogenized 
elasticity tensor is found by volume averaging the constrained elastic properties of the unit 
periodic cell, and these constrained properties, in turn, depend on the elasticity tensor Dijki 
of the homogenized constraining medium. 

The constitutive relations given in equations 126 and 143 are used to update the stress- 
strain history in the constituent phases, whilst equation 150 is used to update the stress-strain 
history in the homogenized self-consistent medium. These relations, which contain Aejj(f), 
A efj(m) and D ljk! , depend on the Eshelby tensor S ljrs for the homogeneous smeared out 
medium, which is defined in equation 123 as 


Sijrs = f 'fj Ujjki (r - r') dV(r') D klrs 
* V 


(170) 


when the field point r lies within the cylindrical volume, V. The “effective” homogeneous 
smeared out medium for a composite with cylindrical fibers will exhibit transverse isotropy if 
the fibers are arranged in hexagonal arrays, and it is shown in Appendix E that the Eshelby 
tensor for a transversely isotropic cylinder, whose X 3 cylindrical axis is normal to the plane 
of transverse isotropy, has the component form. 


£1111 = 

5Dn 1 1 + 71)122 

(171) 

8 D 1 1 1 1 

S-2222 = 

Sim 

1172) 

C — 

3-Cb 122 ~ ^nii 

(173) 

*31122 “ 

8T>im 
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S 2233 

Sum 

S2211 

s 1212 

S 2323 


D 1133 
2 D 1 1 1 1 
S22M 
S 1122 

S 1221 = 


3£>mi — £>1 


122 


8£>im 

S 2 M 2 — -S' 1 3 1 3 = *^1331 = 4 


(174) 

(175) 

(176) 

(177) 

(178) 


with all other S’,,*./ = 0. If the fibers are arranged in tetragonal arrays, the Eshelby tensor 
will exhibit tetragonal symmetry. This case is currently being worked out. 


7 Integration of Self-Consistent Model 

Fourth rank tensors can be written in Voigt notation as matrices and second rank tensors as 
vectors, (c/. Appendix 2 of Mura’s book, [24]). For example, with the notational changes, 


A tT, = A<7,i, 
and 

As 1 = Ae Hi 


A 02 = A(7 22 ! 


As 2 — Ae 22 , 


A<7 3 = ACT;{;{, 


As 3 = As 33 , 


A<t 4 = A cr 23 , 


Ae 4 = 2Ae 23 , 


Acr 5 = A<7 13 , 


Ae 5 — 2Ae 13 , 


Act,; = A<7\2 


Ae§ — 2Aei 2 


Hooke's law for an isotropic elastic medium can be written as 
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A -F 2 fi 
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Ae 2 
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A T 2 fi 
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Ae 3 
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0 
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0 
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Ae 4 

A<t 5 


0 

0 

0 
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0 


Ae 3 

. AVfi ; 


0 

0 

0 

0 

0 

^ . 


. As 6 , 
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For a transversely isotropic medium — such as the smeared out “effective” matrix for hexagonal 
fiber arrays the relationship can be written as 


A<7 1 
Act 2 

A<7 3 

A(7 4 * 

Aar, 

, A<7c, J 


£> mi 

£>1122 

D 1153 

0 

0 

£>1122 

£>1111 

£>1133 

0 

0 

£>1133 

£>1133 

£>3333 

0 

0 

0 

0 

0 

£>1313 

0 

0 

0 

0 

0 

£>1313 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 



' Asj ' 
Ae 2 
Ae 3 

A<f4 

As 5 

. Ae (i . 


(180) 


from which the isotropic results can be recovered by^ taking D I1U = 2/m(1 - is )/ ( 1 - 2 is), 
D.uu - 2/m( 1 - is) /{l - 2is). Du 22 = 2/«V(l - 2 is), Dn33 = WC 1 - 2/'), and D VM * = //. 
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The Eshelby tensor [35] relates the constrained strain increment, A in an inclusion 
which undergoes a transformation or eigenstrain increment, A £* kl , in an infinite medium with 
elasticity tensor in the form 


Aef, ; = SukiAet 


" ij 


In Voigt notation we have 
where the Eshelby matrix takes the form, 


Aef = SijAe* 


5 D 1111 T T? 1 1 22 3T>h 22 ~ -Din i 
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8D l 


in 


8Z?i 


in 


2T>mi 


0 


0 


(181) 


0 


[Sij] = 


3-P 1 122 — -Pi 111 
8 -Dim 

0 


5D U1 1 + -Pi 122 

8-Dmi 

0 


Pll33 

2T>nii 

0 


0 

0 


0 

0 


0 

0 


0 

0 

0 


0 

0 

0 


0 

0 

0 



The integration of the self-consistent model then proceeds as follows. 


0 

0 


/ 3-Pnn — Dim A 

\ 8-Duii / 

(182) 


1. Initialize the starting variables: time t = 0; temperature T = To; overall 
“effective” stress and strain erf = ef = ef = 0 for i = 1 to 6; stress and strain in the 
respective phases <jj(/) = ef(f) = ef(/) = 0, and < 7 j(ra) = ef(ra) = ef(ra) = 0 
for i = 1 to 6; equilibrium stress in the respective phases Sli(f) — Cti(m) = 0 for 
i = 1 to 6; drag stress in the respective phases K(f) — Ko(f) and K(m) = Ko(m). 

2. Compute the overall “effective” elasticity matrix iteratively from the relation 



Y± D f 

v Dik 


s. 


kj 


Yu 

v, 


j_ n ni 

+ u ,k 


+ S U D (D> mJ - + 

[«*; + S u D l l!(DZ,-D mj )]~ 


where is the Kronecker delta matrix, <5*., = 1 for k = j and <5* v = 0 for k ^ j, 
and the Eshelby matrix S, y is given in equation 182. 

3. Start the loading history step. Evaluate the inelastic strain and state vari- 
able increments in the fiber and matrix phases from the unified viscoplastic con- 
stitutive relations. Any unified viscoplastic model may be used. Such relations 
may have the following form: 
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In the fiber the inelastic strain increment is 


Aef (/) = 


( Si (f) - 


(§*,(/) - «,(/)) 

(§*,(/) - «,(/)) 

\ K{f) ) 

K(f) 




where the equilibrium or back stress increment is calculated from the relation 


A «,(/) = elAefif) - 

for / = 1 to 6, and the drag stress increment from 


AK(f) = [qI - e{ ( K(f ) - Ko(fj)\ 

In the matrix a similar set of constitutive relations can be used, so that 


Ae','{m) = 


' Sj(m) — Qi(m) 

K(m) / 


(I «»H ~ ^?(/)) (§s«M - ^M) 


K(m) 


Aa,(m) = g™ Aef (m) - g^^Ae^ (m)Aef(m) Qi(m) 
for i. = 1 to 6, and 


AK(m) = [q™ - q 1 ? (K{m) - /f 0 (m))] \J\ tAe£(m)Ae£(m) 

The quantities n f, n m , g* p , g™, for p = 1 to 4, are material constants associated 
with the unified viscoplastic constitutive relations. 

The deviatoric stress in the fiber phase is defined by 

«*(/) = <b:(/) - | (tfi(/) + <72(/) + < 7 3 (/)) for * = 1,2,3 

and 

Si(/) = cr,;(/) for * = 4, 5, 6. 

Similar relations apply to the matrix phase, viz . , 

s,(m) = <r, (m) - 1 + <r 2 (m) + for i = 1, 2, 3 

and 

Sj(rn) = cq(m) for i = 4,5,6. 

4. Compute the “effective” inelastic and thermal strain increments from the 
relations 

D v (S„ -/„)]" (B„ Aef (/) + C* Aef (m)) 


A£, P = 


and 

where 


a* AT = [D ip (S pq - /„)] _1 (B qk al f + C qk ar) AT 


Bqk- — , . D qj 6jp + S im D mn {D^p D np ^ ( S pk 6pk) 


and 


K 

V, 


0„ = [ft* + S im D m l (£>Z - D np)] 1 ( Sp k - V) 
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5. Evaluate the deviation strain increments from the relations 


and 


Ac q (f) = A<(/) + a*/ AT 
A c q (m) = A ef (m) + a* q m AT 


Ac q = Ae q + a; AT 


6. Evaluate the phase volume averaged total strain increments 
A ef(f) = lh + S ip D-'(D{ J -D, j )}~'{Ae« + 


and 


A ej (m) = 


+ S jq (A c q (f) - Ac,) + S,„D m ' (D{ r - D qr ) Acr(/)} 

+ S.pP„‘ (O” - A»)]"‘ {Ae“ + 

+ Sj q (A c,(m) - Ac,) + (D” - Dr) A<v(m)} 


7. Calculate the stress increments in the fiber and matrix phases from the 
relations 

A Mf) = D{ s (AeJV) - Acj(f)) 

and 

A crj(m) = D™ (AeJ(ra) — Acj(m)) 

8. Compute the overall “effective” stress increment from the relation 

A<r t ° = Da (A ,e° - A if - a* AT) 

9. Update the variables: 


<7i (/, t + At) = 
a i (m, t + At) = 
Uj (/, t + At) = 
f X ( m , t + At) = 
K(/,f + At) = 
(m, t + At) = 
£ i (fyt + At) = 
ef (m, t + At) = 
(/? t + At) = 
ef (m, t + At) = 
of (t + At) = 
£ i (t + At) = 
T(t + At) = 


(/, t) + A(Ji(f) 

(Ti (m, t) + A<7j(m) 

a (/,t) + Aa(/) 

Q* (m,t) + AOj(m) 

K(f,t) + AK(f) 

K ( m,t ) + A K(m) 

£ i (/> *) + Aef (/) 
ef (m, t) + Aef(m) 
£ I (/>t) + Aef (/) 
ef (m, t) + Aef (m) 

(*) + A^ 

e«(t) + Ae« 

T(t) + AT 


10. Start new load step. 
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8 Subvolume Method 


8.1 Approximate Integration of Integral Equations 

The determination of the stress and strain increments throughout the composite material 
requires the solution of the integral equations 

D i]kA4i( r ) = D ?jkA c ki{ r) - 6D m { r) [&e° kl - Ac w (r) - 
- SD ljU ( v)y Y. Y E s«- (0 /// DZ„A<, (*•') e' Ur ~ r,) dV{ r') ( 183 ) 


or 


DZ u Asm = DT jU Ack,(T) - SD m (r) [A4 - Ac„(r>; - 

-«Ai«(r)///£/ tlm „(r-r')Z)“ n „Ae*,(r') dV(t') (184) 

V' 

at each field point r in the unit periodic cell. 

Nemat-Nasser and his colleagues [25,26,27,28] have demonstrated the efficacy of dividing 
the unit cell into a number of subvolumes and assuming that Ae* nn (r') is replaced by 


A4„ M = A E “!, = ~ fJI A4„ (r') dv (r') 


(185) 


corresponding to its average value in the /T h subvolume. 

Let there be N subvolumes in the unit cell, with M subvolumes in the fiber and N — M 
subvolumes in the matrix. Then the preceding integral equations can be written as 

D^kA4i( r) = D? jk Ac k i{r) - 6D ijU (r) [A^ - Ac«(r)' - 


1 ioo ^ r r r \ 

- 6D ijkl (r) — Y £ (0 e‘ (r Y fh- ty dVW) DZnrA< S , 

n v = 0 fl=\ \ J S, J ] 


(186) 


and 


\Vg 


D? ]k AelA) = ^,Ac,Kr)-m (jW (r)[A4-Ac, / (r) 

- W ijk ,{ r) £ £ ff jlhinn, (r - r') dV(r') 

yY; 


D rn Ae* 0 

mnrs rs 


(187) 


/ 

where Vj q denotes the /T h subvolume in the r/" 1 unit periodic cell, and it is assumed that the 
field point r is in the first periodic cell for which q—1. 

These equations can be volume averaged over the a 41 ' subvolume in the unit periodic cell 
to give 


D'r ik AeV> = DVj u Acfc - 6D? jm [Ae‘>, - Acfil - 6D? jU x 


(0 (^fffe^dVir)) £ vJ^fffe-^'dViv^DZnrAe^ (188) 

lr »„o ;*=i \ v e J i J ; 


vW 
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and 


D™ U A4? - D% tI Ac° t ,-6D“ Jk ,{A4,-A4,}- 

6D rn £ £ Jffm r) / // CW, (r - r') rfV'(r') Ae;? 

,=1 ^‘ V " V n v„„ ) 


(189) 


In these equations the deviation strain increments Ac£, are evaluated from the uiiified 
viscoplastic constitutive relation for the cd h subvolume based on the stress value a "(/) or 
crfjirn) in the subvolume, according as the a th subvolume is in the fiber or the matrix phase, 
respectively. The notation 6D^ jkl also denotes the value of D f ljkl — D™ kl or 0 according as the 
a th subvolume is in the fiber or matrix phase, respectively. 

If we use Nemat-Nasser’s notation and write 


Q a (() = ^Jffe' (r dV(r) 


and denote 


Vo 


1 Vc 


(190) 


(191) 


as the volume fraction of the a th subvolume, then the preceding equations may be written as 


N 

£ 

0=1 


iboo 


+ SDf jU E £ £' 9m™ n (0 DZ m f<r «) Q ' 3 (-«) 

n p -0 

- Ac£ 




= £>£ h Ac& - 


(192) 


and 


N 

E 

,i=i 


(r - r') dE(r') 

9-1 V v " 

= ^AcJ, - 6D° m [AeJ, - A<,j 


D m 

mnrs 


Ae*J 


(193) 


where 6 <nj = 1 if a = 0 and S Q & = 0 if a ^ 0, and no sums on a, 0 are intended unless 
explicitly stated. 

Now SD?j kl = 0 if the a th subvolume resides in the matrix. In this case equations 188 and 
189 show that 

Asl‘,' = Ac" for M < a < N (194) 

Thus, only M unknowns (associated with the subvolumes in the fiber) are involved in Ae*f 
and the N — M known quantities (associated with the subvolumes in the matrix) given by 
equation 194 may be taken over to the right hand side of the equations. Equations 192 and 
193 may therefore be written in the form 
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M 

E 

0=1 


.f. rxj 


+ fDfa £ £ E'»...« (O DZ„r‘Q" «) 0 " (-£) 

n p =0 

= A? wAcg “ bDZjki [Ae? ; — Ac?, 


Ac"' 




±oo 

- E E E E' 9«- «) DZnr.fCt (£) <?"(-£) Ac?. 

/i = A/-f-l Tip — 0 


(195) 


and 


A/ 

E 

tf=l 


/ 


A" E ir///^( r )/// A,™ (r - r') dF(r') 

9=1 “ V« V*, 

= D&Acg - SDf su [Ac?, - Ac?,] - 


d: 


A et 


N 


( 


E 6D °m E Ip-/// rfV'(r) /// V Um „ (r - r') rfV(r') 

«=*'+> «-• V " " vi 


D m Ac^ 

Trinrs rs 


(196) 


for a = 1 to M. 

By defining the fourth rank tensor and the second rank tensor ft" in the Fourier 

series representation as 


±OC 


At, = + 6D“ u Y. E E'»«~. (0 D Zmf e 'Q“ «) 0 s (-£) (197) 


and 


6“ = A"«Ac£ - 60g W [A4 - Ac?,) - 

/V ±oo 

- E M>5*< E E E' 9«~» (O DZ„r,fQ a (() Q e (-£ ) Ac?. 


N 

E 

/3=AT+1 n p =0 

or, in the Green’s function representation, as 

l 


(198) 




Eh/// HI u “>™ (r - rfr(r') 

9=1 V '* v;, r*. 


and 




^y«Arjf ; - *£>'} w 


V 

Ae° w - Arfc 

/ 


D 


m 

mnrs 


(199) 


v ^ / \ 
E E h/// d r« // / £W. (r - r') rfV-(r') 
«-■ V " it >>; / 


A’lr.Ac?, (200) 
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the integral equations can be cast in the form. 


M 


£ A%, A £ ;f = Jg for a = 1 to M 

>9=1 


(201) 


In Voigt notation the fourth rank tensor A"^ rs can be written as a matrix A^, and the 
second rank tensors Ae*f and 6" can be written as vectors A ef and 6“, so that 


M 


E = b« for a = 1 to M 

0=1 


( 202 ) 


This represents a system of 6 M linear equations for the unknown values Aef, each matrix 
element af3 of the matrix A consisting of a 6 x 6 submatrix, in the form 


' [An] ••• • 

[Aim] 


{Ae* 1 } ' 


{b 1 } ) 

[A 2 i] • • • 

[A 2 M ] 


{Ae* 2 } 


{b 2 } 

: * ' * [A aj g] * * * : 

i 

{Ae*^} 

> = < 

M 

[Ajvn] • • • 

[Aa/m] 


{Ae* M } j 


, { b "} J 


(203) 


where the submatrix elements are defined as 




1 

£ P 

A <*0 
/1 12 

A 1 a/3 

"913 

A a 0 

-9 14 

A a 0 

A a/3 

9]6 

ta0 

/i 2 l 

AO.0 

^22 

ACC0 

^23 

4 a/3 

4 a/3 
-9^25 

4 a/3 
-926 

AOC0 

■^31 

A<X0 

/1 32 

AOC0 
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-'934 

A a / 3 
9-35 

A a/3 
-936 

\CX0 

^41 

4 a/3 
•'9-42 

A <*0 
^*43 

4 a/3 

^44 

A a / 3 
945 

A a/3 

946 

4 a/3 
^51 

Aa0 

1 952 

4 a/3 
^53 

4 a/3 
^54 

A q 0 
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A a/3 
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4 a/3 
-^61 

A a/3 
/1 62 

A <*0 
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965 

4 a/3 
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(204) 


and the corresponding column vectors as 


{Ac*"} = 


Aef ] 




6f 1 

Aef 




6? 

Aef 

Aef 

► 

and 

M = • 

*>3 

1 

Aef 




# 

. Aef , 




bg J 


(205) 


This system can be solved by standard Gaussian elimination. However, if M subvolumes 
are included in the fiber, these equations represent a 6A/ system of equations, whose solution 
may pose storage problems on the computer. 
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8.2 Solution of Integral Equations by Iteration 

An alternative is to use equations 188 and 189 in an iterative fashion. As a first guess the 
integral terms in 188 and 189 can be neglected and we obtain 

D? jk Ael? = DT jk A4i - SD? jU [A 4 - A cfi] for a = 1 to M 

corresponding to the subvolumes in the fiber, and 

D% kl Aet? = D% kl Ac a kl for a = M + 1 to N 


corresponding to the subvolumes in the matrix. These relations can then be substituted 
into the integral terms to yield an improved “Rayleigh-Born” approximation to D ™ kl for 
a = 1 to N . This process can then be repeated until D™ kl Ae* k f converges to within a user 
specified tolerance. In essence we solve the equations 


V?m } a+1 = D mA<, - [AeJ, - Ac?,' - 

M ±oo 

- E SD !ju E E EW (0 DZ,r«fQ" «) a 6 (-€) {as;?}, 

i^=l ?!/>=() 

N ±oo 

- E E E S ’ (0 c,/<r (£) H) Acf s 

(3=M + 1 rip— 0 


/*=A/ + 1 


q= 1 \ 


(206) 


A*, { A£ ir} A+1 = A"«Ac£ - [a £ » - Acj] - 

- E W Oi)U E I A Iff dV( r) ff fu Umn (r - r') dV(r') ) AU { Ae* } - 

0-. I" \ ” 

- E dDf ju f i {^JjfdV(v)f]Ju Um Ar-v')dV(v')\DZ nr A4, (207) 


until the (A + l) th iterate differs insignificantly from the A lh iterate. 

In the solution of the composite problem, two constituent phases, namely the fiber and 
matrix phases, have been considered. For composites with a third chemically degrading phase 
separating the fiber from the matrix, the preceding solutions may be modified by assuming 
that, in the summations from (3 — 1 to M, some of the subvolumes, say from f3 = L to 
A/, pertain to the degraded material. It will then be necessary to postulate a viscoplastic 
constitutive relation for this chemically degrading phase. 

The total strain increment in the a u ‘ subvolume in the unit cell is then obtained by 
averaging equations 57 and 80 over the ct th subvolume to give, 

M ±oo 

A el? = A 4, + Y. EEE'»« (O Ar (O QC (-0 r>z,rAt4 + 

1 n r — 0 

+ E EEE'sh (C)/‘ i Q“«)0''(-€)£’™„Ac» (208) 

■ M -t I 0 
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or 


A4“ = A £ », + t t ( £ ///' dV(r)j f f u tlm „ (r - r') dV(r') 

TV oo ( 

+ £ £ 


\ 


D m Ae*'* + 


V, 


Pq 


£= A /+ 1<?=1 


vIII dv(l) III Uk,m " (r “ r<) c!l,(r ' ) 

“ Va V* / 


/ 


£) m Ac' 3 

mnrs rs 


(209) 


for a = 1 to AT. 

The constitutive relation, required to update the stress and state variables in each sub- 
volume, is then given for the a th subvolume as the average of equation 35 in the form, 

A<7“ = £>”« (A e{? - As;,") (210) 

If we assume that N = 2, with one subvolume in the fiber and the other in the matrix, 
then the theory is similar to the self-consistent model in which the strain increments in 
the constituent phases are assumed to be spatially constant and equal to their respective 
constituent volume averages. However, the interaction effects of the nearest neighboring 
cells are fully accounted for since geometric periodicity is assumed in the integral equation 
formulation and the material outside the unit cell has not been smeared into an “effective” 
uniform material. 

In both the Fourier series and the Green’s function formulations integrals of the form 


JJJ e iK r dV (r) (211) 

v;, 

need to be evaluated over the subvolume, V a . These Laue interference integrals [39] can be 
evaluated exactly if each subvolume consists of a circular or oblong cylinder. In the case of a 
circular cylindrical fiber, each subvolume within the fiber would consist of an infinite cylinder 
with a cross-section in the shape of an element of area in cylindrical coordinates, comprised 
of two circular arcs with constant radii, rj and r 2 , and two radial segments along the lines 
of constant 6\ and d 2 - An attempt will be made to evaluate equation 211 for this type of 
cross-section. If this proves too unwieldy, the subvolumes within the cylindrical fiber can 
be taken to be cylinders themselves, with the cylindrical fiber represented as a “bundle of 
sticks” . We assume that the actual fiber is comprised of subvolumes of the correct shape, but 
we make an approximation in performing the volume integration over a circular cylindrical 
subvolume. 


9 Concluding Remarks 

This document is the first annual report on NASA Grant NAG3-882. Much of the work on 
which this report is based exists only as a melange in the literature and we have therefore 
attempted to write the report in enough mathematical detail that it can be worked through 
without reference to the literature. In the second year we shall work out the required integrals 
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in the formulations and program the methods in FORTRAN subroutines suitable for inclu- 
sion in nonlinear finite element programs. In the third year we will determine the material 
constants for various composite materials and provide a comparison of the present theory 
with finite element and experimental results. 

Our aim is to produce an end product which can be used in nonlinear finite element 
and boundary element programs for analyzing the structural behavior of composite materials 
under thermomechanical loading conditions at elevated temperature. 

The viscoplastic behavior of periodic composites is analyzed by means of implicit integral 
equations. These integral equations arise when the problem of determining the stress-strain 
variation throughout a unit periodic cell in the periodic composite is solved by a Fourier 
series or Green’s function approach. In this report we show that the Fourier series and 
Green’s function approaches are mathematically equivalent by means of the Poisson sum 
formula. By applying simplifying assumptions the integral equations can be solved in an 
approximate fashion and used in structural analysis programs to analyze the overall behavior 
of the composite. When the strain-temperature history at the “damage critical” location 
has been determined from the structural analysis, this can be used to “drive” the “exact” 
integral equations to determine the stress-strain history variation throughout a unit periodic 
cell located at the critical location. 

The unit cell in the periodic structure can be formulated to analyze fibrous, laminated 
and particulate composites. By retaining the effects due to the application of displacements 
and tractions at the surface of the composite it is also possible to analyze the behavior of thin 
walled composite sections such as are found in turbine engine combustor liners and blades. 
When this is done the integral equations which must be solved are basically those which are 
used in boundary element programs. In the constitutive subroutine which we plan to embed 
in the nonlinear finite element program to analyze the overall macroscopic behavior of the 
composite, we effectively have a boundary element equation (specialized for the case of a 
periodic composite) which we solve in an approximate fashion for the stress at the Gaussian 
integration point when the boundary displacement on the element is prescribed by the finite 
element program. 

When the effects of damage are included in the constitutive formulations it will be possible 
to embed the subroutine in an optimization program such as ADS in order to determine 
optimum composite configurations. 
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Appendix A 

Properties of the Green’s Function 

Consider a point force /*, ( r') acting at the point r' in an infinite medium with elasticity tensor 
Djjiei. From the definition of the Green’s function the displacement at the field point r due 
to the point force f k (r') at r' is 


«i(r) = Gij (r - r') fj (r') 
so that the infinitesimal strain at r is 

1 ( dG tj (r - r') dGrnj (r - r') ’ 


^m(r) = - 


dx. 


+ 


dxi 


fj (r') 


and the associated stress is 


or 


^kp ( I* ) Dkpim 2 


Gkpi?) 

1 ( dGij (r - r') dG mj (r - r') 


dx 1t 


+ 


dxi 


fj (r') 


(A.l) 


(A. 2) 


(A. 3) 


(A. 4) 


Since the elasticity tensor D kpim is symmetric with respect to the indices i and m , the last 
relation can be written as 


_ r> dG V ( r “ r ') t tJ\ 
^kpi^) — Dkpim p. fj (r ) 

OXm 

For static equilibrium, we must have 

fk (r') = - jf n„(r)a pk (r) dS( r) 


(A. 5) 


(A. 6) 


where S denotes any closed surface in the infinite medium with an outward unit normal n p (r) 
which surrounds the point of application of the point force f k (r'). An application of Gauss’ 
divergence theorem gives 


M = -/ff^W--fff°~ 


d 2 G v (r - r') 


T IJ 

9Xj yidx p 


Writing 


fk (!•') = fj (r') fff ShjS (r - rj dV( r) 


fi(vjdV(r) (A.7) 


(A. 8) 


then gives the equilibrium requirement that 

d 2 Gj j (r - r') 

I \ ^kpim 

V ^ 


Miff 


dx m Ox,, 


+ S kj S (r - r ; ) [ dV (r) = 0 


(A. 9) 
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(A. 10) 


Since f 3 (r') and V are arbitrary, then 

f) 2 G (r - r'\ 

D kpi m + M ( r - r') = 0 

is the differential relation satisfied by the Green’s tensor function. When multiplied by fj, 
this is just Navier’s equation of elasticity with the displacement Uj(r) = G l} (r — r') fj (r') and 
the body force set equal to 6 k jfj (r') 6 (r — r'). 

Rearranging the indices, this differential relation can be expressed as 

= 0 (A - n) 

The solution to the differential equation can be found by applying Fourier integral tech- 
niques. On multiplying the differential relation by e lhqXq dx\ dx 2 dx 3 , i.e., by e lK r dV( r), and 
integrating over all space, we obtain 


Dij k i J j J ^ ^ k Q~ e lhqXq dx i dx 2 dx 3 + S ip J j j 6 (xx) 6 (x 2 ) 6 ( x 3 ) e lKqXq dx x dx 2 dx 3 = 0 


-d 2 G kp ( r) ^ 

V 3 

^V/ 

(A. 12) 

From the sifting properties of the Dirac delta function the last integral is unity, so that 



Da 
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ijkl 


+ Di 
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{ dGkpjr) ' 
dx,\ \ dxj 

( dG kp (r) \ \ 


+ D 


ijk2 


( dGkpjr) 


d 


dx 2 \ dxj 


+ 


e lKqXq dx i dx 2 dx 3 + 6 ip = 0 


'i?7c3 Vv 1 o 

dx 3 \ dxj j j 
Integration by parts severally with respect to xi,x 2 ,x$ then gives 


(A.13) 



dG kp (r) 
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Dij k l£ 

x e i{KxXx+Kzxi) dXi dxs + 
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J x 1 = — OO 


Dtjfc3e 
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1 K^dG kp {r) 


JK,X? dG kp {r) 
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Dij k2 & 
e i(K ' xl+K2x2) dxldx2 


dxi 

3 J X2~ — OO 
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X3 = -00 


}- 


-///{ 

T bin — 0 


n — KVK ~ , ifC 1 D -iv- 1 n 1 iK a x : 

Dijki — — Um + Dij k 2 — — iK 2 + Uijkn — 7 ^- — e 


q dxi dx 2 dx-i + 
(A.14) 


The surface integrals are zero since dG kp (r)/dxj vanishes at the infinite lower and upper 
limits of integration, so that one integration by parts yields the result, 

- D ijkl JJf iK '~i~-- elKr dV ( r ) + S iP = 0 (A. 15) 

-OO J 
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A second integration by parts yields 


OO 


+ D m HI i 2 K l KjG kp (r)e iKr dV( r) + 6 ip = 0 

— OO 

(A.16) 

or 

Dij k iKiKjG kp (T£) Si P 

(A. 17) 

where 

OO 


G kp { K) = I 1 y^Gfc P (r)e iK,r dV (r) 

— OO 

(A.18) 

is the Fourier integral transform of G kp (r). By writing £ as a unit vector in the direction of 

the wave vector K, we have 

r _ Kj _ Kj 

1 <]K,K, K 

(A.19) 

in which K = K q K q is the magnitude of the K vector. Then 


D,, u ^^-K 2 G tr (K) = S ip 

(A. 20) 

or 


K 2 Dij k i(iCjG kp (K) = 6i P 

(A.21) 

The Christoffel stiffness tensor M (c/. [33]) is defined by the relation 


Mik{ C) = D ijklClCj 

(A. 22) 

so that 


K 2 M ik (C)G kp (K ) = 6 ip 

(A. 23) 

Premultiplying both sides by the tensor K~ 2 M -1 gives 


Gij(K) = K~ 2 M^(C) 

(A.24) 

The Fourier inverse of equation A. 18 gives 


OO 

G,j( r) = (2n)- 3 JIJe- Kr G ij (K)d 3 K 

— OO 

(A. 25) 

where d' ! K = d,K t dK? dKj , so that we finally obtain the Green’s function 
integral form 

in the Fourier 

r fri [[[ ( ^K Mij '(O -/K.r 

G ' Al) JJJvn) 3 A ' 

— OO 

(A. 26) 

with 


«.-■«) = (Awoa) -1 

(A. 27) 
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This representation of the Green’s function yields explicit results for isotropic and transversely 
isotropic materials (c/. Mura’s book, [24]). For cubic and general anisotropy the Fourier 
integral representation must be used. 

Often, we are concerned with volume integrals of the Green’s function and its derivatives 
with respect to r, such as C//t< mn (r). It is then advantageous to use the Fourier integral 
representation even for isotropic and transversely isotropic materials. The advantage is gained 
by reversing the order of the wave vector and volume integrations, whereby many of the 
integrations can be carried out explicitly. 

Sir William Thomson (Lord Kelvin) obtained an explicit form for the Green’s function of 
an isotropic elastic material in 1848. As an example we may deduce the Kelvin result for the 
Green’s function of an isotropic material from the Fourier integral relation. For an isotropic 
material 

Dijki T "b &u&jk ) (A. 28) 

and so M t k(C) = Dij^QQ has the form 

M ik ( C) = (ACiCfc T ^SikClCl + lK\ iCk) = ((A + flKiCk + fttiik) (A. 29) 

since 

C/C/ = Ci+C 2 2 + C3 2 = (^) +(^) +(^) =1 (A. 30) 

The inverse tensor is given by the relation 

M lk \0 = ((A + riC * + /■**)'' = S f - (A.31) 

which is easily verified by showing that 


Mij'Mjk - 6 ik 


(A. 32) 


From the preceding relations 


MZ l M jh 


' 6ij A + 


H p(A + 2/i) 


OCj J ( ft)CjCk + 


x + |S (A + /x) 2 (A + //)// 

'SiU- + o xGsfc , r, \SlSk — Oik 




n(\ + 2/r) p(A + 2/j,) 


as required. 

The Green’s function is therefore obtained in the form 


G„(D = (*>)■* Jff-£ e -iKr d 3 K— (2 -)-/// (A + 2/() A , 

— -oo 

From a table of Fourier transforms (c/. [36]) we find that 


A + ^ ^^ p -,K.r d 3 K 


(A. 33) 


(A. 34) 


(A. 35) 
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which may be differentiated with respect to x t and Xj to give 


d 2 r 


1 d > K 


K‘ 


dxidxj JII 

-OO 

By contracting the i and j indices we obtain 


The Green’s function may therefore be written as 


GiA r) = 


1 / 8 a , d 2 r 


V rrr 2 

-7T 


A + fi 


-7T 


<9 2 r 


or 


(27r) 3 y // dx q dx q /x(A + 2/i) dxidxj 
A + /i 


Gii(r) ~~ 87r//r { 2 ^ ** A + 2^ r 2 ' ) } 


where the relations 


d 2 r 


dx q dx q r 


and 


A d" 2|ti 

d 2 r 


8ij X{Xj 


dxidxj r r 3 


obtained by differentiating r = ^/x q x q , have been used. 


(A. 36) 

(A. 37) 

(A. 38) 
(A. 39) 

(A. 40) 
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Appendix B 

Relationship Between Fourier Series 
and Green’s Function Approaches 

In the composite material the total strain increment Ae^(r) is periodic in r and is defined 
by the relationship 

Aeju(r) = Ae° kl + Ae kl (r ) (B.l) 

where Ae kl is the strain increment applied to the composite’s boundary which is equal to 
the volume average of Ae kl (r) over the unit periodic cell, and Ae k[ (r) is the deviation or 
perturbation from the average value due to the presence of the fibers. 

From equations 84 and 85 the perturbed strain increment is given in the Fourier series 
and Green’s function approaches by the equivalent relations, 

Ae«(r) = £'»«« (0 /// D£M„ M dV( r') (B.2) 

C n P=° V, 


l(r) = III Uut ’ (r " r '> D Zr,^r, M dV(r') 


We now show that these equations are equivalent and that the Green’s function relation is 
the Poisson sum transformation of the Fourier series relation. 

From the definition of g k i mn ( C) i n equation 54 we may write 


9kiij ( 0 = \ (M ik l (C) CjQ + MaHOCiCk) 


(BA) 


9kUj ( Ci i C2 ) C3 ) — \ (M ik l (( 1, C2, C 3 )CiCi + M u 1 (Ci,C 2 ,C 3 ) 0 a) (B.5) 


where 


c- = - 
c ' e 


We may therefore write 


2nn-> \ 2 


(no sum on i ) for i — 1,2, 3. (B.6) 


wj( 0 = 9kiij( Ci (nun 2 ,n 3 ) ,C 2 (n u n 2 ,n 3 ) ,C 3 (n u n 2 ,n 3 )) = f klij (n u n 2 , n 3 ) (B.7) 


The perturbation strain increment can then be written in the form 


boo ±00 ±00 


AeiJr) 


L]L 2 L 3 


LI? fklij (ni, n 2 , n 3 ) / ff DV] r , Ae* rtl (r') 


ti \— 0 11-2— 0 713=0 
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or <us 


where 


±OC i 00 .too 


Ae w (r) = ----- E E E' hki(ni,n 2 ,n 3 ) 

"1-^2 "3 m=0 n-2= 0 713=0 


(B.9) 


hti{n,,n 2 ,n 3 ) = fkuj (n u n 2 ,n 3 ) JJJ D," ,Ae;, (r') x 

/27T71W x 27m 2 f . x 27T713 . , v\ 

xe ‘( ( X1 - X i) + l 2 (^-x 2 )+ Ls (*3-* 3 )j da .; ( B . 10 ) 

By the Poisson sum formula (c/. Morse and Feshbach’s “Theoretical Physics”, [37]) 

±oo ±oo ±oo ±00 ±oo ±oo 


E E E'Mni,n 2 ,n 3 ) = E E E w/// d3 Ke^'^^x 

. _n « o n — — 0 m 3 =0 • / _‘oo‘ / 

K 2 L 2 K 3 L 3 


77i=0 7l2~0 713=0 


X few 


2tt ’ 2tt ’ 2t r 


0 


(B.ll) 


where the sum over the integers ni,n 2 ,n 3 is replaced by the sum over the integers mi,m 2 ,m. 3 
in the Fourier integrals. The sum over m, includes the case where m\ = m 2 = m 3 = 0. 

We now have the alternative sum, 


±co ioo ±oo 


-j XtV 

Ae«(r) = - — — E E E' ^ (ni,n 2 ,n 3 ) 

12 3 7 / j () r? 2=0 713=0 


±00 ±00 ±00 


_i_ ^v; _l_ vV a, ^ /» /* n Jo 1/ 

= S s E 


+m 2 K2L2+m 3 K3L3) 


fkl 


771] =0 777 2 = 0 7713=0 Eqq 

777 . . 

X 


K x Li K 2 L 2 K 3 L 3 


(2tt) 3 ~ JKUJ \ 2tt ’ 2tt ’ 2?r 

JJJ D™ rs Ae* s (r') dx; ds^da^ (B.12) 


or 


**> - £££///» * • X 


mi =0 777 2=0 777-3=0 ' 


x JJJ D% r „As* ra (r') e _ *( Al dx) dxf 2 dx^ 

lv 

(B.13) 

Due to the geometric periodicity of the unit cell we may write 

A e*. s (r') = Ae* s (x'^ x'.^ x^) - As,*,, - mxL u x' 2 - m 2 L 2 ,x 3 - m 3 L 3 ) 


and 

dx, dx 2 dx a — d (x\ — m\L\) d (x' 2 — m 2 L 2 ) d (x 3 — m 3 L 3 ) 
so that by making the change of variable 

(x) - r/i]Ti, x' 2 — rn 2 L 2 ,x' a — m 3 L 3 ) = (x" , x 2 , x") = r" 


(B.14) 

(B.15) 

(B.1C) 
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we obtain 


A£ «< r > “ jJlS^4^^^) e,iK ' n+K ~' 


X 


±oo ioo ±oo r r r 

* £ £ £ / / / 

7711 =0 7712=0 7713=0 f/ / ^ \ 

Vc (mi, m 2 , m 3 ) 


A7r«^ £ rs ( r ") e~ iK r " dV (r") (B.17) 


where the volume integration extends over the volume K (m\,m 2 ,mz) of the unit cell whose 
center is at the point ra 2 L 2 , m^L^). Since mi,m 2 , 777,3 range over all integer values, 

the summation of the volume integrals extends to all the cells in the periodic lattice, i.e ., 
it extends over the entire volume, V, of the composite medium. The expression for A £ki{r) 
thus takes the form 


A£ «< r) = 


(2ir)= 


K\L\ K 2 L 2 K 3 L 3 \ ;K.r 
27r ’ 27r ’ 27 r 


///A>, Ae;, (r") e~ IK •" dV( r") 


(B. 18) 


By interchanging the order of the volume and wave vector integrals and noting that r" 
can be replaced by r' since it is a dummy integration variable, we obtain 

A £a( r) = Sjjtm ///H im (^. ^r) £,K 


Introducing ( 


V -00 

/Cl Lj /C 2 L 2 K 3 L 3 
27T ’ 2n 5 27t 


) in place of (ni,n 2 ,ri 3 ) in the expression for 


fkiij (ni , n 2 , n 3 ) = 9kiij ( Ci (n x , n 2 , n 3 ) , C 2 K , n 2> n 3 ) , Cs («i , n 2 , n 3 ) ) (B.20) 


then gives 


9klij Cl 


K1L1 K2L2 K3L3 

2tt ’ 27T ’ 27r 


> C 2 


K X L X K 2 L 2 K 3 L 3 

27 r ’ 27 t ’ 27r 


1 C3 


KxLi K 2 L 2 K 3 L 3 

27T 27T 27T 


1 ' 


(B.21) 


with 


Ki = Kj 

K s/KrK, 


(B.22) 


and the perturbed strain increment takes the form 


OO 

Ae ui r)=JIJdV(r')JJJ 


d 3 K 1 / M^(0 
(27 t) 3 2 \ K 2 


KjKi + 


jgno 

X 2 


p ' K ■ ( r - r ' ) s Ae^jr') 
(B.23) 
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But. from Appendix A, 


* <* - «o - — - ///s^ — * (“> 


(2tt) 3 K 2 ~ JJJ (2tt) 

— CXJ — OO 

since Ga- (r — r') — Gik (r' — r), and therefore 

S 2 G j4 ( r-r')_ fffd’KMi'K) 


- -Ill i 




dxjdxi JJJ (27 r) 3 K 2 

J -30 

Inserting the last relation into the expression for Asjt;(r) then shows that 

,, 1 /Wi* (r - r') t 3 2 G« (r - r') 


Ae H (r) = -///^(r')- 
1/ 


+ 


D™ Ae* fr') 

^iirs ^rs V / 


dxjdxi dxjdxk 

From the definition of the tensor Ukimn (r — r') in equation 83, we see that 

Ae«(r) = jfj U mi (r - r') £>” , Ae'„ (r') dV(r') 


(B.25) 


(B.26) 


(B.27) 


which is the result obtained with the Green’s function approach. 

The Fourier series expression for the perturbation strain increment is thus identical to the 
Green’s function expression and the two are linked via the Poisson sum formula. 
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Appendix C 


Poisson Sum Formula 

In Fig. 7 the function f(x ) = [(a; — is shown on the unit cell extending 

to x = L. The corresponding function defined on the n th unit cell to the right 
f(x + nL ) and the periodic function q(x), which is comprised of the functions 
defined on all of the unit cells extending from x = — oo to X = + 00 , is given by 


from x = 0 
is given by 
f(x + nL) 


OO 

Q(x) = J2 f( x + nL ) 

71 = — OO 


(C.l) 


Each function, /(x + nL), is defined only over the corresponding n th periodic cell and is taken 
to be zero outside of the cell. Each function can therefore be represented as a Fourier integral 
and the periodic function q(x) can be written as a sum of Fourier integrals, 

OO OO 

q(x) = 5Z f( x + n L) = 5Z Fourier integral of f(x + nL) (C.2) 

n — — oo n — — oo 


By setting x = 0 in both summations we obtain the Poisson sum formula. This method is 
outlined at the end of this Appendix. 

The Poisson sum formula can also be derived by expanding the periodic function q(x) into 
a Fourier expansion and showing that the Fourier integral sum, when x = 0, is the sum of 
the coefficients in the Fourier series expansion. 

Since q(x) is a periodic function of period L, it may be expanded into a Fourier series in 
the form 


°° s i2irmx 

q{x) = J 2 a m e ~ L 


(C.3) 


m~— oo 


where 


am ~ L l 

The object is to show that the Fourier series 


L i27rmx / 

e l q (x') dx' 


(C.4) 



i2nmx 

e~ l 



i2nmx' 

e l q ( x ') dx' 


(C.5) 


represents a sum of Fourier integrals. This is easily accomplished by introducing the expres- 
sion 

OO 

Q ( x> ) = J2 f( x> + nL ) 

n=— oo 

into the Fourier expansion and changing the integration variable by means of the relation 


y = x' + nL 


In this way we obtain 


OO 

q(x) = f(x + nL) 

71 — — OO 



j 27T771X 

L 


E 


71 — — OO 


/ 


y=(n-l)L 


i2n7n(y—iiL) 

e r f(y) dy 


(C.6) 


(C.7) 


52 



The exponential function exp(i2Trmy / L) is a periodic function with period L, so that 


i'2nm(y—nL) i2irmy 

e l = e ~ 


(C.8) 


If we set x = 0 and note that the sum over the integration limits is equivalent to summing 
over the entire axis of x from x = — oo to x = +oo, we obtain 


^ w 1 poo i2nmy 

9(0)= E f(nL)= E 7 / e 1 f(v) d v 


Putting L = 1 gives 


OO OO -oo 

9(0)= £ /(n)= £ 7 e‘ 2nmy f{y)dy 

, „ J —00 


(C.9) 


(C.1U) 


and by changing the integration variable to K = 2ixy/L, we obtain 

(KL' 


OO OO T rOQ / U I \ 

z f{n)= z^L eimKLf (-^) dK 


(C.ll) 


In three dimensions this result takes the form 

±00 ±00 ±00 ±00 ±00 ±00 


E E E/K.«2.«3) = E E E 

ii=0n2 — 0 713=0 mj=0m2=0 rri3=0 v —00 

A1L1 K2L2 K3L3 




27T 27T ’ 27T 


0 


(C.12) 


which is the form used in Appendix B. The cubic function defined here for illustration purposes 
has the property that the constant ao in the Fourier expansion is zero, since f 0 L f(x) dx = 0. 
This term may therefore be omitted from the summation on the left and the summation signs 
primed to denote the omission of the term with rii = n-i = n 3 = 0 . 

It is now possible to show that the Poisson sum formula follows from the Fourier integral 
sum in equation C.2. 

We have the Fourier integral sum representation 


OO OG -oo 

E /( m ) = E / f(v)Hy - m ) d v 

m= -oo m= — oo* 00 

where the Dirac delta function is given by the Fourier integral 

✓ 

6(y — m) = — r e i(y ~' n)z dz 

271 " ./ — 00 

Then the Fourier integral of f(m) is 

1 POO POO 

f(m) = — e~ imz dz / e izy f(y) dy 

Z7T J — 00 J —oo 


(C.13) 


(C.14) 


(C.15) 
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Putting 2 = 27 ra gives 


/ OO roo 

e -i2*ma da ^™Vf(y)dy (C.16) 

-oo J -oo 

The Fourier series expansion which culminates in equation C.10 shows that 

/(<*)= fe'-’/H-l! (C.17) 

J — oo 

so that equation C.16 becomes 

f(m)= r e~ i2 * ma f(a)da (C.18) 

J — OO 

We may therefore write 

OO -OO 

Y /M = Y f( m ) 

m =— oo m— oo 

This is the Poisson sum formula in 


OO 00-00 

= E /(-"»)= E / e‘ 2 ’~> /(<*)<*« (C. 19 ) 

m= — oo m=— oo ^ 00 

equation C.10, from which equations C.ll and C.12 follow. 
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Appendix D 


Integral Equation for Displacement Increment 
In Neighborhood of Free Surface 


The static equilibrium equation for a medium with elasticity tensor D™ kl 
equation 37 as 


D 


m 

ijkl 



_d_ 

dxj 




is obtained from 


(D.l) 


in which A^(r) = Aeg,(r) + Ae w (r). 

In this equilibrium relation we will not assume that the strain increment Aejy(r) applied 
at the surface of the composite is constant and will take it as a spatial variable. If we set 


A/, (r) = ~ {d”«a £ ;,w} 

(D.2) 

and note that 


. r / \ 1 ( d ( Au k(r)) d(Auf(r))\ 

A4(r) = 2( ' dx, ’ + \x t ’) 

(D.3) 

the equilibrium equation may be written in the form 


<1 

II 

V £ 

<1 Sl 

v — " H 

Q 

(0.4) 


where the symmetry of D™ kl with respect to the indices k and l has been used. On denoting 
the operator Tn by the relationship 


the equilibrium equation is 


Now consider the integral 


r - D" 1 92 

~ v '’ a 9x~ir t 


FaAuf = Afi 


(D.5) 

(D-6) 


1(4,, V>) = (r') r„4>j (r') d,V( r') (D.7) 

V' 

for any two field variables <p,(r) and ipi( r). These field variables may be tensors of any rank. 
For example, if </> and ip were second rank tensors, then ip) would be a second rank 
tensor integral 

Ipqifr tp) = Jff <P P i (r') Fijipjq (r ; ) dV (r') (D.8) 
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From the definition of the operator T t j we have 


!(*,*) = f/J M*>) os* £ 7 


(D.9) 


and since 
d 


7 I <X M D m 


dx\ 


dipj (r')\ _ dfc (r') dj) 3 (F) 


AA D m .muji + <£. ('F') 

, ^ipqj f)„' ^ V 1 v 1 > u 


ivqj dx’ 


d ( dipj (r')’ 


-p \ --'g 

the integral becomes 


ax; ~ tw ax; ' wj ax; V ax; 


(D. 10) 


'<**> ^ ///j£j (* < r ') (D H) 

The first integral can be transformed into a surface integral via Gauss’ divergence theorem, 
so that 

m,i>) = // (D. 12) 

s n v P q 

By interchanging the arguments <f> and t/> it is evident that 

W, 0) = // ", ( r ) 0i M O™ 0S(r') - lll d -^DT ml d -^dV(r') (D.13) 

*5 <7 y P Q 

Now the elasticity tensor jD^ is symmetric with respect to its indices, so that the interchanges 

ip -> pi, qj -*■ jq, ip +-> ( qj or jq), qj «-> (ip or pi) (D-14) 

leave the elasticity tensor unaltered. This shows that the volume integrals in /(</>, i/>) and 
/( t/i. 0) are identical, so that Green’s identity ([38], page 434) can be written as 

/«>,</>) -/(0,0) = //», M (* M - * M c^^r) (D is) 


or 


III (tSiAi - *<?**,) dV = II n p U,D - ADtJM ds 

\7 Q \ Q Q / 


(D.16) 


Now choose the field variables <fi and ip as a vector and second rank tensor in the forms 

(f>i (r') = A uf (r') - A u° (r') = A a, (r') (D.17) 


and 


P-,j (r') = Gij (r - r') 


(D.18) 


where 1 (p, (r') is the perturbation displacement increment Aw, (r') and Gj, (r — r') is the Green's 
tensor function satisfying the differential equation (cf. Appendix A), 


F,.iG, k (r - r') + 6 ik b (r - r') = 0 


(D.19) 
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The integral relation then becomes 


fff {(At if (r') - A (r')) TijG jk (r - r') - 

V 

- G* (r - r') Tij (AuJ (r') - Au» (r')) } dV(r’) 
= jjn r (r') I (Auf (r') - A«» (r')) £>: 


aG jfc (r - r') 


m 

ipqj dx' 


G lk (r - r') DT mj 


'd(Au J (r')) _ d(Au° (r ; )) 

dx’ 


dx' 


dS( r') 


(D.20) 


'9 ~~9 

Replacing TijGjk (r — r') by —8 ik 6 (r — r') in the first term of the volume integral gives 
Att^r) - Au° k (r) + JJJ G ik (r - r') T {j (A uj (r') - A (r')) dV( r') 


- JJn p (r') 


V 


Gik (r - r') D? nj 


( 8(AmJ (r')) a(Au°(r')) ’ 

\ dx’ t 9x' q 


(Au[ (r') - Au° (rj) D™/ 


9Gjk (r - r') 


dx' 


dS( r') 


(D.21) 


From the definition of the operator Tj we obtain the relation 

XijAuJ(r) = A/,(r) = {Z^A^r)} 

and by inserting this result into the integral equation and noting that 


(D.22) 


D 


d(Au$ (r')) 


tpv dx' 


= D m - 
2 


w q - \ 9 

we find that the total displacement increment is 

d 


( d(Au° (r')) ^ d(Au° (rQ) ^ _ 

V dx 9 


dx' 


= D m Ae°- (r') 
ipqj qj v x ) 


(D.23) 


Aul(r) = Aul(r)-JJfG ik (r-r')£ F {D™ a (Ae;A^-Ae° rs (r'))}dV(r’) + 

v j 

+ Jf n v (O { Gik (r - r') D™ mj (AeJ (r') - Ae° qj (r')) - 

C v 


(Auf(r')-Au»(r'))Aw : 


( r - r ') 


dx' 


dS( r') 


(D.24) 


In the volume integral an integration by parts with Gauss’ divergence theorem using the 
relationship 

dG ik (r - r') _ dG ik (r - r'_) ^ 25 ^ 


dx] 


dxi 
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gives 


r - r') ~ {A>, (*C (r') - AC (r'))} dV'(r') 


= Jjj~{G tt (r-r')D? Jr ,( A £ ;,(r')-A £ * (r'))} dV(r') - 

V' 3 

- r ' ) flS,. (ACM - ACM) dK(r') 

V’’ 3 

= JJ nj (r') G',* (r - r') (Ae^ (r') - Ae° s (r')) dS{ r') + 


+ 


in 


dG ik (r - r') 

dxj 


D” , (AC (r') - A4 (r')) dV(r') 


(D.26) 


This result may now be substituted into the integral equation to produce 

AuJ-(r) - A<(r) - (AC M - AC (r')) dV(r') + 

V 3 

+ JJ n, (r’) ( G* (r - r') D”, (AC (r') - AC M) + 

5 ^ 

+ (A u r(r')-A„»(r'))AC aG ' t(r_r ' ) ' 


dx r 


dS{r') 


(D.27) 


In the first two terms of the surface integral we observe from equation 35 that 

n, («•') D? jr , (AC (r-) - AC M) = n. ,A<r« (r') = At, (r') (D.28) 

represents the incremental surface traction on the surface of the composite. Equation D.27 
represents the well known Somigliana identity ([38], page 93) for the displacement increment. 
In the case where the composite is assumed to be of infinite extent the surface integrals in 
the preceding integral equation vanish, and if As®, (r') is assumed to be spatially constant, 
the total displacement increment is given by the relationship, 

A*) = A«»(r) - /// a °‘ ~ - D?ir, AC M dV( r') (D.29) 

V ‘ 3 

which corresponds to equation 79 and is the form used in the main report. However, equa- 
tion D.27 must be used when the surface is not infinitely removed and if As®, (r') is not 
assumed to be constant. 

In a finite element context it will normally be assumed that the fibers are very small in 
comparison with the dimensions of the finite element. At the Gaussian integration point in 
the finite element it is then permissible to neglect the contribution from the surface integrals 
since the surface of the finite element is assumed to be many periodic cells away at “infinity’' . 
In some situations, however, this may not be a valid assumption. Some turbine blades and 
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turbine engine combustor liners are fabricated from thin sections in which the central passages 
are hollow to allow cooling air to pass through the component. In the thin cross sections of 
such components the surface integrals must be retained in the constitutive formulation. 

Suppose, for example, that the total displacement increment at the node points of a finite 
element are given. From these nodal values and a knowledge of the element’s displacement 
interpolation functions it is then possible to compute the total displacement increment Au°(r) 
on the surface of the element and the total strain increment Ae° s (r) at any point. Since 
A uf (r') = A a? (r 7 ) on the surface of the finite element, the last term in the integral equation 
vanishes and the total displacement increment is determined from 


Auf(r) 


A«2(r) - jjj dG>k( J x r) A7r. K M - A4, (r')) dV( r') + 

V j 

+ JJnj (r') G it (r - r') (A£ (r') - A e' r , (r')) dS( r') (D.30) 

s 


in which the terms in the surface integral represent the contribution to the total displacement 
increment due to the incremental traction, 

At, (r') = rij (r') D"j rs (A£ (r') - A £ ;„ (r')) (D.31) 

on the surface of the element. This surface traction is needed to maintain the displacement 
increment equality A uf (r 7 ) = Auf (r 7 ), which is imposed at the element’s surface. 

By differentiating Au^(r) with respect to and xi and taking half the sum, the total 
strain increment is subject to the integral equation 


A4(r) = Ae» (r) + JJJ U mi (r - r') D%„ (AeJ. (r') - A£ (r')) dV( r') + 

+ fjnj(r') 


1 / dGik (r - r ; ) dGji (r - r’) 

2 l dxi dxk 


D™ , K, (r') - A £ ;, (r')) rfS(r') 

(D.32) 


in which 

OJlAeli (<•') = A”V.,Ac« (r') - SD ijtl (r') [a £ J (r') - Ac t , (r')] (D.33) 

and this integral equation should be used for thin sections of composite material where sur- 
face effects are important. This implicit integral equation is similar to that for the infinite 
medium but contains a correction term for the surface effects in the last integral. This surface 
integral will become less important — due to the derivatives of the Green’s function ---when 
the integration points r' are far removed from the field point r and it vanishes for an infinite 
medium. 

In the preceding development it was assumed that the displacement increment A uf (r') was 
known, by interpolation with the element displacement polynomials, from the nodal values. 
This forces the incremental surface traction At, (r 7 ) to adopt a periodic distribution in order 
to maintain A uf (r') = A n® (r 7 ) on the surface of the element. We could, alternatively, assume 
that the surface traction increment is zero on the free surface of the element, in which case 
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the total displacement increment Au[(r) will exhibit a periodic variation on the the surface 
and the surface takes on the appearance of a frilled structure. 

If we therefore assume that the finite element is thin (see Fig. 8); that the surfaces are 
free of surface traction; and that the surfaces at the ends of the finite element are sufficiently 
far removed from the Gaussian integration point, the first term in the surface integral in the 
integral equation is zero and in lieu of equation D.30 the relationship for the total displacement 
increment now takes the form, 


A<(r) 


A„2(r) - /// aGi g" - Pgr. (A<v, (r') - Ae r “, (r')) dV( r') + 

V j 

+ jjn, (r') (Auf (r') - Au° (r')) O'", — ^ ~ ^ dS( r') (D.34) 


The solution to this integral equation gives a periodic total displacement increment, Au[(r), 
which, on the surface of the composite, will exhibit frilling. 

It is clear that during the finite element analysis frilling will not occur in the element. The 
interpolation functions normally used in isoparametric elements are linear and quadratic, and 
cannot adopt the required periodic behavior. However, the stiffness of the finite element — as 
computed at the Gaussian integration points with the composite constitutive model — will 
reflect that the fact that the constitutive properties are computed as though the element were 
free to take on a frilled appearance. When the “damage critical” strain-temperature history is 
used to determine the stress-strain history variation throughout the unit periodic cell outside 
of the finite element program, the preceding integral equation will allow the frilled appearance 
of the composite to be calculated. 


60 



Appendix E 


Evaluation of the Eshelby Tensor 


The Eshelby tensor S ip i m is defined by the relation 

= "5 { &Ss/// G " (r - r ° dV{t,) + 6^kM G -’ <r " r ° dV{l,) ) 3 ^ (E1) 

or as 

Sipim = JJJ U iprs (r - r') dV(r') D rstm (E.2) 

v 

where the field point r lies within the volume, V, and where the volume extends over an 
infinite cylinder of radius a in a medium with elasticity tensor D ljk i. Although the Green’s 
function for transversely isotropic materials is known [24], it is more convenient to work with 
the Fourier integral representation of the Green’s function as given in Appendix A. 
Introduction of the Fourier integral representation, 


OO 

G “ < r - r ' J = /// 

— OO 


d 3 K 

(2tt) 3 K 2 


e -i K.(r-r') 


(E.3) 


where = Ki/K = K x j ^JK q K q , into one of the volume integrals in the definition of S ip i ri 
gives, on reversing the order of the volume and wave vector integrations, 


Lkyij — 


d 2 


dx k dx g 


///G«(r-r')c(V(r') 


(E.4) 


or 


_ d 2 f f f d J K 

Lkgij ~d^ g JJJ ( 2 ^ 


d 3 K A/,/(C) -t’K.r fffjK.r 


(27r) 3 K 2 


JJf^ iK r ' dV(r') (E.5) 

v‘ 


The Laue interference integral [39] extends over the cylindrical volume and can be written as 

1 = JJJ e iK v dV{ r') = JJJ e '(^\+^x' 2 +K 3 x' 3 ) dx ^ dx t^ dx ^ ( E . 6 ) 


Let x\ gcosO, a:, = psin#. Then in cylindrical coordinates 

r oo ra r 2n 


/ oo ra rliK 

/ / e <A'3* 3c »(A'.tfcostf+A' 2ff 8 "'Vdx'zQdQdO 

-oo Jo JO 


Since 


/ OO 

e iK * x * dx' A = 2 tt8 (K s ) 

-OO 

where 8 (Ks) is the Dirac delta function, the integral takes the form 


/ = 2ir8(K :i ) f a f 2 ’ eHKie<*»o+K 3 g*in &) gd Qdd 
Jo Jo 


(E.7) 

(E.S) 
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Let q = 


, dg = 


dq 


sJKl + Kl y/Kl + Kl 


. Then, if K = yflcf+K f, 


If we now set 


( <?/CiCos0 g/C^sinfl \ 

V^+^2 + V^+^ 2 ) qdqdO 

Jki + k% j » y ° 


(E.IO) 


cos 


= sin#', then 


y/Kf+K$ sJK[+kI 

I = -MS /°* / 2?r g dq dQ 

sjK\ + K\ J o Jo 


(E.ll) 


Since the integration extends over a whole circumference, it is immaterial where the origin 
of 0 is placed. The integral may therefore be written as 


/ = 


2tt< 5 (tf 3 ) f lK 


& iq cos 6 dQ 


or as 


^HkIJo qdq Jo ' 

•4=^=/ qdq2TTJ 0 (q) 
y/Kf + K\ Jo 
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J, 
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where Jo and Ji denote the usual Bessel functions of order zero and one. 
The integral Lkgij can therefore be written as 


(E.12) 


(E.13) 
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(E.16) 


If A- = 3 or q = 3, the Dirac delta function 6 (/\ 3 ) gives zero values for the integral. Hence, 
the non-zero values of are given by A: = 1, 2 and p = 1,2. 
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Invoking the sifting properties of the Dirac delta function, viz., 


/ OO 

/ (K u K 2 , K 3 ) 6 (I<,) dK 3 = f (K u K 2 , 0) 

-OO 


then gives 


= -if] dK ' iK ‘ M « «■>&•& = °> {Hfk) 

where the unit vector C is now defined by the relations 

Kx , K 2 


( K \ X \+ K2X - 2 ) q 


0 = 




C2 — 




C3 — 0 
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Jx (ay/Ki + K?) 
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(E.18) 


(E.19) 


If we put 


Cl = -TT 


K x 


COS#, C2 




Ko 


K y/Kf+Kl ’ K y]K* + ATI 

and set Xx = r cos (j), x 2 = r sin (j), then in cylindrical coordinates, 


sin# 


(E.20) 


Lkgij ~ 


2tt 


1 K dK d0 M u (Ci , C2) CkCg e ~ iKr «« 9 -+ ) (E.21 ) 


The integration with respect to 9 extends over a complete circumference, so that 

1 r 2lT 


Lkgij 


(E.22) 


J p Ztt r 00 
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Since S>: p ; m is real, the real part of the preceding integral involving the integration over K is 

r 00 

TV = / a cos(A> cos #) J\(aK ) dK (E.23) 

Jo 

Setting z = rcos#, and noting that cos (Kz) = cos(— Kz), we need be concerned only with 
positive values of 2 . Now if the field point r lies within the cylindrical volume, then 0 < z < a. 
But, from Gradshteyn and Ryzhik [40], 


N= r a cos ( K z)J\ (aK) dK = 

Jo 

If tv — sin" 1 (z/a), then 


o cos {sin x (z/a)} 


N = 


acost/i 

V^T 2 


\/a 2 — 


cos ^ cos xl> 


for 0 < z < a 


(E.24) 
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Thus, 

U m = «i.6)f*C > de ( E - 26 > 

independent of position r in the cylinder as expected from Eshelby’s result. In this integral we 

have Ci = cos 0 , C2 = sin < 9 , C3 = 0, A/," 1 ((1,(2) = (CmAmjnCn) and k and g are restricted 
to the values 1 and 2. The Eshelby tensor may now be written as 


Siplm = -^Dkjlm | ^ Mpj 1 (Cl) C 2 ) CiC h d0 + M tj 1 (Cl) C 2 ) CpCk d# j (E. 27 ) 

When C:i = 0 the Christoffel stiffness tensor for a transversely isotropic material, M tj , and its 
inverse, M t ]\ (which applies to the homogenized medium of a composite with fibers arranged 
in hexagonal arrays) have the component forms 
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The Eshelby tensor can now be determined by integration in the form, 
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b I I 33 — S?2Xl 

S-221] = <S'n22 
«Sl2]2 = ‘S'l221 = 


8A lil — A 122 


8A111 

S2323 = S2332 = •S'lsai — ‘S'1313 = = 3 


(E,14) 

(E.45) 

(E.46) 

(E.47) 


The Eshelby tensor for tetragonal materials -which applies to the homogenized medium 
of a composite with a square array of fibers — is currently being worked out. 

The results for an infinite isotropic cylinder may be recovered by taking 


film = 2p(l — u)/(l — 2v), D u 22 = 2nv/(l -2v), and Am = 2fxu/(l - 2v) (E.48) 


where p is the Lame shear modulus and v is Poisson’s ratio. For an infinite isotropic cylinder 
the Eshelby tensor reduces to 


5 - \v 
5,111 8(1 — u) 

(E.49) 

S2222 = S\m 

(E.50) 

C 4 / 2-1 

5,122 8(1 — v) 

(E.51) 

q _ 

1 2233 2(1 - u) 

(E.52) 

Si 133 — S2233 

(E.53) 

S 2211 = S1122 

(E.54) 

s - s - 3 ' 4u 

5,212 - 51221 - 8(1 - u) 

(E.55) 

S2323 — ‘S'l 313 = “Smi = S2332 = \ 

(E.56) 


The Eshelby tensor for both isotropic and transversely isotropic materials can also be 
deduced from equations 17.27, 17.30 and 17.31 of Mura’s book, [24], by setting g = 0 in his 
notation. 
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Appendix F 

Proof that ^(x - y) = £%/( y - x) 

From the definition of U t j k i(x — y) we have 

TT , \ 1 (9 2 G ik (x-y) , d 2 G jk (x-y)' 

Uljkt(x y) 2 [ dxjdx t + dxidxt , 

But 0G *[ X - y) = - aGit l x ~ y) , so that 


dxi 


dyi 


d 2 G ik (x - y) d 2 G ik (x - y) 


dxjdxi 


dyjdyt 


The operator can therefore be written as 


TT , , 1 (d 2 G ik {x-y) d 2 G jk {x-y)' 

tWx - y) - - j ( ' -g-y^T + 9y,9y,"\ 


But G ik (x - y) = G ik ( y - x), so that 


U, jk i(x - y) = -- 


1 / d 2 G ik (y - x) d 2 G jk { y - x) 
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2 \ ify/diA ' %% 

Z/»jw(x - y) = U ijk i( y - x) 


as required. 
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Appendix G 

Differentiation of Singular Integrals 

In the text and Appendices we have taken derivatives of the volume integrals and written, 
for example, 


hq — 


d_ 

dx q 

III 


III 


dG ik (r - r') 


dx. 


D™ s Ae* rs (r') dV( r') 


(r-F) 

dx q dXj 


D™ s Ae* rs (F) dV( r') 


(G.l) 


If the integration volume V contains the field point r the integrand dGik (r — F) /dxj is 
singular at the point r' = r, and the above operation in which the derivative is taken inside 
the integral must be treated with caution, as pointed out by Bui, [41] and Born and Wolf, 
[42]. We should, in fact, isolate a small spherical volume, D, about the singular point F = r 
and evaluate the integral according to Bui’s procedure, viz., 


h. 
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-III 
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dx q dxj 
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0 yr) 
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where we have used the fact that, if the spherical volume D about the point r is small enough, 
the strain increment can be considered constant and taken to have the value at the center of 
the sphere, A£* s (r). The integral may therefore be written as 


Ikq — 


III 

V-D 


d 2 G lk (r - r') 


dx q dxj 


DTirAK. X) dV( r') - 


II n ,X) 


,, dG ik (r - r') 
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dS(r') D”' r ,Ae;. s ( r) 


(G.3) 


The first volume integral is evaluated in the principal value sense as D — > 0. 

Rather than using the preceding operations outlined by Bui, we may treat Gjj (r — r') 
as a Fourier integral. The preceding operations are not then necessary and the derivative 
can be taken inside the integral. That is, equation G.l is valid when the Fourier integral 
representation of the Green’s function is used. 
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To demonstrate the validity of equation G.l, consider the singular integral used by Bui. 
He considers the derivative of the integral 

f 1 dt 

F(x) = / = log 

J - 1 t — x 

where — 1 < x < 1. Since the integral is known, its derivative is simply found as 


1 — x 


1 4 - x 


(G.4) 


dF 

dx 


x — 1 X + 1 


(G.5) 


Notice that the integrand is singular at the point t = x. Bui demonstrates that in order to 
take the derivative of the integral we must write it in its principal value sense, 


Fix) = lim 

v ’ f^O 


r r dt 

J t — x J t — x 


(G.6) 


t=-'l 




and the derivative dF/d,x must be evaluated by noting that both limits and the integrand are 
functions of x. Using Leibnitz’s rule for differentiating an integral whose limits depend on x 
gives 


dF 

dx 



(G.7) 


or 


dF 

dx 


1 


(G.8) 


x — 1 X + 1 

To avoid the converted terms which arise from differentiating an integral whose limits 
depend on x, consider representing the integrand as a Fourier integral. We have, from Grad- 
shteyn and Rhyzik [43], the Fourier integral representation, 


= ( G ' 9 > 

The singular integral F(x) can then be written as 

F{x) = L dt wJ-J 

= TsTJ rwmeV'ilC (G.10) 
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If we now differentiate with respect to x in the normal manner we obtain 
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A comparison of this integral with equation G.9 shows that this Fourier integral has the 
inverse relation, 


dF _ __L 1_ 

da; x — 1 x + 1 


(G.12) 


which is the correct result. 

Thus, by expanding the integrand of a singular integral as a Fourier integral, reversing the 
integrals, taking the normal derivative, and inverting the resulting Fourier integral, we obtain 
the correct derivative of the singular integral. It is then clear that if the Green’s function is 
represented in Fourier integral form the procedure of Bui is not required. In fact, the Eshelby 
tensor in Appendix E is obtained by taking the derivative of the Fourier integral, and the 
correct result is obtained. 


69 


Appendix H 

Origin of Self-Consistency 

Many researchers in the mechanics literature suggest that the self-consistent method has its 
origins in the present century. It would appear that the method is, however, very old and 
has its origins in the last century. In the Lorentz-Lorenz theory ([42], page 87 and [45]) of 
1880 the electric dipole moment p in a dielectric is related to the electric field E' by the 
constitutive relation p = evE', where a is the polarizability. The polarizability a is related to 
the refractive index n and the number of molecules per unit volume, N. If E is the mean or 
volume averaged field applied to the dielectric the actual field at any point is given by 

47T N 

E' = E+-^p (H.l) 

o 

where 47rA r p/3 denotes the pert urbation or deviation from the average electric field. As shown 
on page 85 of reference [42] this value is estimated by smearing the effects of the molecules 
outside a spherical volume enclosing the point at which the field is observed. An analogous 
formula for statical fields had been derived even earlier by Clausius in 1879 and Mossotti in 
1850. 

Twersky [44] observes: 

In the biography of John William Strutt (third Baron Rayleigh) by his son Robert 
John (the fourth baron), the son quotes the father on the verse that faces the ini- 
tial contents page of the first four of Lord Rayleigh’s six volumes of Scientific 
Papers: “When I was bringing out my Scientific Papers I proposed a motto from 
the Psalms, ‘The works of the Lord are Great, sought out of all them that have 
pleasure therein’. The Secretary to the Press suggested with many apologies that 
the reader might suppose that I was the Lord.” The Secretary need not have been 
so apologetic. The second verse of Psalm 111 should have been augmented with 
the next three lines: “His work is honourable and glorious, and his righteousness 
endureth forever. He hath made his wonderful works to be remembered.” Depart- 
ing from King James' translation, we may read in the Hebrew of the last verse 
of this psalm the most important of all the Rayleigh principles of mathematical 
physics that the wise beginning of work in this field is to assume that the prob- 
lem had been considered by Rayleigh and to study his works: “The beginning of 
wisdom is reverence for the Lord; very good sense have all who do so.” 

Rayleigh [45] tackled the problem in his paper “On The Influence of Obstacles Arranged 
in Rectangular Order Upon the Properties of a Medium” and was probably the first person to 
define when the self-consistent method, viz. the Lorentz-Lorenz formula, could be expected 
to break down. At the end of his paper he states: 

The general conclusion as regards the optical application is that, even if we may 
neglect dispersion, we must not expect such formulae as (the Lorentz-Lorenz equa- 
tion) to be more than approximately correct in the case of dense fluid and solid 
bodies. 
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FIGURE 2. - UNIT PERIODIC CELL FOR LAMINATED 
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FIGURE I). - FLOWCHART OF FINITE ELEMENT IMPLEMENTATION. 
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FIGURE 5. - PERIODIC UNIT CELL IN HEXAGONAL FIBER ARRAY SURROUNDED 
BY NEAREST NEIGHBOR CELLS. 



SMEARED OUT 
"EFFECTIVE" MEDIUM-; 


/ 



FIGURE 6. - PERIODIC CELL IS REPLACED BY A SURROUNDING MATRIX OF 
"EFFECTIVE" RADIUS b AND THE NEIGHBOR CELLS ARE REPLACED BY A 
SMEARED OUT "EFFECTIVE" MEDIUM WHOSE AVERAGE CONSTITUTIVE PROP- 
ERTIES ARE THE VOLUME AVERAGE OF THE CONSTITUTIVE PROPERTIES 
IN THE FIBER AND MATRIX, WHEN CONSTRAINED BY THE "EFFECTIVE" 
MEDIUM. 



UNIT CELL OF LENGTH L, FROM x = 0 TO X = L. THE PERIODIC 
FUNCTION q(X) EXTENDING FROM - oo TO oo IS OBTAINED BY SUMMING 
THE FUNCTIONS ON THE UNIT CELLS IN THE FORM 


oo 

q(x) = J^fCx + nL) * 
n = -oo 



POINT 


* THESE SURFACES ARE FAR ENOUGH AWAY FROM THE 
GAUSSIAN INTEGRATION POINTS TO IGNORE THE 
SURFACE TRACTION CONTRIBUTION TO THE TOTAL 
DISPLACEMENT INCREMENT Au! (D 


FIGURE 7. - REPRESENTATION OF A PERIODIC FUNCTION AS A SUM OVER FIGURE 8. - REPRESENTATION OF A FINITE ELEMENT FOR THIN COM- 

UNIT CELLS. POSITE SECTIONS WHEN SURFACE IS FREE OF APPLIED TRACTION. 
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